Genome-wide association study in individuals of European and African ancestry and multi-trait analysis of opioid use disorder identifies 19 independent genome-wide significant risk loci

Despite the large toll of opioid use disorder (OUD), genome-wide association studies (GWAS) of OUD to date have yielded few susceptibility loci. We performed a large-scale GWAS of OUD in individuals of European (EUR) and African (AFR) ancestry, optimizing genetic informativeness by performing MTAG (Multi-trait analysis of GWAS) with genetically correlated substance use disorders (SUDs). Meta-analysis included seven cohorts: the Million Veteran Program, Psychiatric Genomics Consortium, iPSYCH, FinnGen, Partners Biobank, BioVU, and Yale-Penn 3, resulting in a total N = 639,063 (Ncases = 20,686;Neffective = 77,026) across ancestries. OUD cases were defined as having a lifetime OUD diagnosis, and controls as anyone not known to meet OUD criteria. We estimated SNP-heritability (h2SNP) and genetic correlations (rg). Based on genetic correlation, we performed MTAG on OUD, alcohol use disorder (AUD), and cannabis use disorder (CanUD). A leave-one-out polygenic risk score (PRS) analysis was performed to compare OUD and OUD-MTAG PRS as predictors of OUD case status in Yale-Penn 3. The EUR meta-analysis identified three genome-wide significant (GWS; p ≤ 5 × 10−8) lead SNPs—one at FURIN (rs11372849; p = 9.54 × 10−10) and two OPRM1 variants (rs1799971, p = 4.92 × 10−09; rs79704991, p = 1.11 × 10−08; r2 = 0.02). Rs1799971 (p = 4.91 × 10−08) and another OPRM1 variant (rs9478500; p = 1.95 × 10−08; r2 = 0.03) were identified in the cross-ancestry meta-analysis. Estimated h2SNP was 12.75%, with strong rg with CanUD (rg = 0.82; p = 1.14 × 10−47) and AUD (rg = 0.77; p = 6.36 × 10−78). The OUD-MTAG resulted in a GWAS Nequivalent = 128,748 and 18 independent GWS loci, some mapping to genes or gene regions that have previously been associated with psychiatric or addiction phenotypes. The OUD-MTAG PRS accounted for 3.81% of OUD variance (beta = 0.61;s.e. = 0.066; p = 2.00 × 10−16) compared to 2.41% (beta = 0.45; s.e. = 0.058; p = 2.90 × 10−13) explained by the OUD PRS. The current study identified OUD variant associations at OPRM1, single variant associations with FURIN, and 18 GWS associations in the OUD-MTAG. The genetic architecture of OUD is likely influenced by both OUD-specific loci and loci shared across SUDs.


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].
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 nonopioid 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 (N case = 5435 N effective = 20,032), a total EUR N = 554,186 (N case = 15,251; N effective = 56,994), and an overall N = 639,063 (N case =20,686; N effective = 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. Cohortspecific summaries of AFR and EUR OUD subjects are presented in Table 1 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 metaanalyses were performed. Measures of cross-sample heterogeneity (Cochran's Q, I 2 ) 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).
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 × 10 −08 ) OPRM1 variants in low LD (r 2 < 0.1).

SNP-heritability and Linkage-Disequilibrium (LD) Score Regression
EUR OUD GWAS summary statistics were used to estimate SNP-heritability (h 2 SNP ), and to characterize OUD genetic correlations (r g ) 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 traitspecific 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 × 10 −05 (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.
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 −08 ), an intronic variant. Rs1799971 was also GWS in the cross-ancestry metaanalysis (p = 4.91 × 10 −08 ), and is not in strong LD with rs9478500 (EUR r 2 = 0.03; AFR r 2 = 0.002; ALL r 2 = 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 crossancestry GWAS (p = 6.00 × 10 −08 ).
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 (r g = 1.02; p = 2.38 x 10 −214 ) [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.14 x 10 −47 ) [12] and AUD (r g = 0.77; p = 6.36 x 10 −78 ) [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.15 x 10 −03 ) [10].
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 r g = 0.98 (p = 1.22 x 10 −77 ) [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 2 = 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 2 = 0.02). Two OPRM1 variants were also GWS in the cross-ancestry OUD GWAS (rs1799971 and rs9478500; r 2 = 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 2 = 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 2 SNP) of 12.75% (Z = 11.28) compared to the previous largest OUD GWAS (h 2 SNP = 11.30%; Z = 6.27) [7]. However, the comparison between these two studies is not direct: the largest previous GWAS [7] used opioid- Table 3. Genome-wide significant (p ≤ 5.00 × 10 −08 ) Lead SNP associations in OUD-MTAG analysis (of 441 total GWS SNPs). 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 nonsignificance 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].