GWAS of cannabis use disorder (CUD) summary statistics on unrelated genotyped individuals of European ancestry (N _cases  = 14 080, N _controls  = 343 726) were obtained from the Psychiatric Genomics Consortium (PGC; https://pgc.unc.edu). These summary statistics are the combination of samples from the PCG Substance Use Disorders working group, iPSYCH and deCODE. The iPSYCH cohort consists of individuals born in Denmark between 1981 and 2005. CUD cases were defined using ICD10 codes (F12.1‐12.2). ¹⁴ Controls were individuals who did not have ICD10 codes related to CUD. DeCODE cases were drawn from the largest addiction treatment centre in Iceland, the SAA‐National Centre of Addiction Medicine. CUD diagnoses in this treatment cohort were made by clinicians using the Diagnostic and Statistical Manual of Mental Disorders (DSM) system (DSM‐IIIR, DSM‐IV and DSM‐5 criteria). ¹⁵ The cases in the PGC Substance Use Disorders working groups all met DSM‐IV diagnostic criteria.

Gene‐based association analysis for CUD GWAS summary statistics ⁶ was performed using MAGMA version 1.09 (https://ctg.cncr.nl/software/magma). MAGMA maps SNPs to genes by aggregating common variants (MAF > 0.01) at the gene level, which increases discovery power by using the linear combination of SNP‐wise P‐values as test statistic, reducing the burden of multiple testing correction seen in univariate GWAS. We used the 1000 genomes phase 3 European reference panel population for LD estimation and mapped variants to the NCBI hg19 genome assembly, which contained 18 297 autosomal protein‐coding genes with SNPs mapped within the defined coordinates of 5 kb upstream and 1.5 kb downstream. Any genes from the major histocompatibility complex (MHC, chr6:28477797–33448354) on chromosome 6 were removed. Multiple testing correction was done using the Bonferroni method, where we divided the alpha level by the number of comparisons and set P < 2.7 × 10⁻⁶ as the P‐value required for significance.

Given the association signal may be therapeutically actionable and to further support the MAGMA findings that PDE4B as the most significant gene in this locus, we subjected the PDE4B region to further analysis to evaluate whether it is a causal association. Firstly, we leveraged probabilistic finemapping to prioritise putatively causal genetic variation in this region. Specifically, we used a conventional Bayesian approach that was applied to all variants in the CUD GWAS within 3 MB of the defined PDE4B genic boundaries. Asymptotic Bayes' factors (ABF) for each SNP were approximated using Wakefield's method, assuming a prior variance of 0.2, ² as outlined extensively elsewhere. ¹⁶ ABF were summed to define credible sets given that Bayes' factors are proportional to the posterior probability (PP) for causality for each variant. In other words, to define a 95% credible set of variants in this region, which contains a causal variant with 95% probability, variant‐wise PP were summed in descending order until 0.95 is exceeded. This method assumes a single causal variant so that we did not have to account for LD between variants, which has been demonstrated to be susceptible to false positives in finemapping studies with a prior of multiple causal variants that use references external to the GWAS or not otherwise very well matched at a population level. ¹⁷

We then wished to further investigate the phenotypic relevance for the finemapped CUD association signal in PDE4B. For the variant with the highest PP, we performed a phenome‐wide association study (pheWAS) using IEU open GWAS project database (https://gwas.mrcieu.ac.uk/phewas/) and the FinnGen release 7 resource (https://r7.finngen.fi/about). The concept underlying this is to assess other traits to which this variant is linked.

We then attempted to infer using expression quantitative trait loci (eQTL) data whether upregulation or down‐regulation of PDE4B was associated with liability to CUD and whether CUD and PDE4B expression displayed statistical colocalisation. PDE4B eQTLs were sourced from the multi‐tissue eQTL catalogue resource (https://fivex.sph.umich.edu/). ¹⁸ We used cortical eQTL data assembled by the MetaBrain consortium (N = 2970) for colocalisation analyses via the coloc method. ¹⁹ The coloc approach infers the PP of five competing hypotheses (H) for a given region: H ₀ = the region is associated with neither trait, H ₁ = the region is associated with trait one, H ₂ = the region is associated with trait two, H ₃ = the region is associated with both traits but with a different underlying causal variant, and H ₄ = both traits share a causal variant. This was implemented using default priors via version 4 of the coloc R package. We also performed a sensitivity analysis where we varied the SNP priors between (1 × 10⁻⁸ and 1 × 10⁻⁴), as implemented by the sensitivity function in the coloc package.

To further refine our understanding of CUD associated genes that may present opportunities for repurposing, we performed a transcriptome‐wide association study (TWAS) and a proteome‐wide association study (PWAS) via the FUSION framework. ²⁰ We achieved this through leveraging genetically imputed models of mRNA and protein expression. For TWAS, we used tissue samples of whole blood (GTEx v7) and post‐mortem brain (GTEx v7, PsychENCODE), ²⁰ , ²¹ whereas protein expression weights of post‐mortem brain and whole blood were yielded from Religious Orders Study and Memory and Aging Project (ROS/MAP) ²² and NHLBI's Atherosclerosis Risk in Communities (ARIC), ²³ respectively. As this method integrates SNP effects from the model of genetically predicted expression with the effects of the same SNPs on CUD, TWAS/PWAS Z‐scores, after accounting for linkage disequilibrium (LD), can be a conceptualised measure of genetic covariance between mRNA or protein expression of the gene and the GWAS trait of interest. As a result, the sign of the TWAS/PWAS Z score is informative as to which direction of genetically predicted mRNA or protein expression is associated with increased odds of CUD. To ensure we captured only the most confidently associated genes that could constitute drug repurposing candidates, we utilised a conservative method for multiple‐testing correction whereby the Bonferroni methodology was implemented to divide the alpha level (0.05) by the total number of significantly cis‐heritable models of genetically regulated expression (GReX) tested from any brain tissue considered or whole blood. Finally, we implemented probabilistic finemapping of TWAS Z scores using FOCUS v0.6.10 to refine potential causal genes for which predicted expression is associated with CUD, as outlined extensively elsewhere. ²⁴ , ²⁵ Using the default prior and prior variance, we estimated marginal posterior inclusion probabilities (PIP) for membership of 90% credible set for each gene in the region in and around the implicated hyaluronidase gene cluster given its potential therapeutic relevance.

LD‐score regression analysis (LDSR; v1.0.1) (https://github.com/bulik/ldsc) ²⁶ was used to estimate genetic correlation between CUD and GWAS on 50 biochemical traits from the UK Biobank (UKBB) as outlined previously by our group for other psychiatric disorders. ¹³ As the mode of action of many existing drugs involves modulating biochemical traits, for example, lipids and blood glucose, shared biology that may be indexed by genetic correlation could be informative as drug repurposing opportunities. In LDSR, the genetic covariance is estimated by regressing SNP‐wise χ², the product of the marginal SNP effects from both traits, on its LD score. The SNP heritability estimates for both traits are used to normalise the genetic covariance to obtain genetic correlation (r _g), which can be accurately estimated in the presence of any sample overlap only affects the LDSR intercept and not the slope. Bonferroni multiple testing correction was used for the 50 biochemical traits tested to define a significant r_g. The CUD and biochemical trait GWAS summary statistics were cleaned to ensure they contain HapMap3 SNPs outside the MHC with minor allele frequency >0.05 for consistency. To evaluate evidence of a causal relationship, we used the latent causal variable (LCV) on all genetically correlated CUD and biochemical trait pairs as demonstrated elsewhere and in our previous work. ¹³ , ²⁷ To estimate partial genetic causality, the LCV framework leverages the bivariate genome‐wide distribution of marginal SNP effects on both CUD and each of the biochemical traits and outputs the posterior mean genetic causality proportion metric (GCP), with GCP > 0 implying partial genetic causality of trait one on trait two, and vice versa. To calculate the GCP metric, the LCV model utilises the genome‐wide SNP–trait association Z scores for two traits and the mixed fourth moments (cokurtosis) of the respective distributions to assess whether there is evidence for a causal effect of one trait on the other. To guard against false positives, partial genetic causality was defined using the recommended threshold of a significantly non‐zero|GCP| > 0.6. ²⁷

Using aggregated gene‐level association (MAGMA), we observed two genes that were associated with CUD—PDE4B (P = 2.09 × 10⁻⁶) and FOXP2 (P = 9.30 × 10⁻⁷) after Bonferroni correction (P < 2.7 × 10⁻⁶; Table S1). These genes have been previously reported in Johnson et al. ⁶ Further analysis using the Drug Gene Interaction Database (DGIdb) on the MAGMA significant genes for CUD revealed the highest interaction score for PDE4B was with dyphylline, a vasodilator and bronchodilator agent. Dyphylline is a phosphodiesterase 4 inhibitor clinically used for the prevention of bronchial asthma or other respiratory diseases. ²⁸ It has similar pharmacological actions and safety profiles as other xanthine derivatives, such as caffeine and theobromine. ²⁹ According to DrugBank, the mechanism of action of dyphylline in humans is a cAMP‐specific 3′,5′‐cyclic phosphodiesterase 4A, 4B, 4C and 4D inhibitor, and an adenosine receptor A1 and A2a antagonist. ³⁰ , ³¹ , ³² , ³³ No drug interactions for FOXP2 were observed in DrugBank or in DGIdb. Interestingly, the PDE3 inhibitor, dipyridamole, also targets PDE4B and has been investigated therapeutically for bipolar disorders. ³⁴

Using probabilistic finemapping, we prioritised one variant in the 95% credible set with a posterior probability greater than 40% (rs1392816, PP = 0.534), out of a total of 34 variants within the credible set. This variant was an intronic variant in the PDE4B gene but was not associated with PDE4B mRNA expression in any of the eQTL catalogue studies or in the MetaBrain cortical eQTL dataset. A pheWAS in the IEU GWAS database of this variant revealed that it was associated with mostly inflammatory [e.g. leukocyte counts, C‐reactive protein (CRP)] and anthropometric (e.g. fat mass, waist‐to‐hip ratio) measures using a conventional phenome‐wide significance threshold (P < 1 × 10⁻⁵) (Tables S2 and S3), However, it also demonstrated an association with smoking behaviours. In the FinnGen GWAS database, which only considers binary phenotypes collected from linked hospital inpatient records, it was associated with alcohol use disorder and other related substance abuse phenotypes, supporting its phenotypic relevance for CUD. Statistical colocalisation between PDE4B expression in the cortex and CUD supported the third colocalisation hypothesis that this region is associated with both PDE4B expression and CUD, but both traits are influenced by a different causal variant (PP _H3  = 0.93). This conclusion remained consistent with different priors. Indeed, there was no strong evidence to support the CUD‐associated variants mapped to PDE4B as acting as eQTLs, suggesting that the effect of genetic variation in this region is mediated through another modality or it influences expression in a tissue or cell type not considered.

We utilised TWAS and PWAS methods to integrate genomic information into functionally relevant units that map to genes or proteins and their expression. Next, we assessed the statistical associations between CUD and predicted gene expression. Although TWAS/PWAS associations do not necessarily imply a causal relationship, this approach can be implemented to identify candidate genes located at loci with an inferred mechanistic underpinning through expression. ³⁵ In multi‐brain region and blood TWAS of CUD, significantly down‐regulated genetically proxied expression of the hyaluronidase gene HYAL3 was associated with the disorder in the whole blood, as well as several brain regions including the amygdala (Z = −4.294, P = 1.75 × 10⁻⁵), anterior cingulate cortex (Z = −5.103, P = 3.34 × 10⁻⁷), hippocampus (Z = −5.103, P = 3.34 × 10⁻⁷), hypothalamus (Z = −5.185, P = 1.05 × 10⁻⁷), nucleus accumbens (Z = −5.552, P = 2.82 × 10⁻⁸) and putamen (Z = −5.433, P = 5.51 × 10⁻⁸). The TWAS also revealed significantly up‐regulated expression of NAT6, also known as HYAL1, in the cortex (Figure 1A).

NAT6 is a hyaluronidase that is proximally located on chromosome 3, ³⁶ similar to HYAL3. However, HYAL1 is considered to be a more canonically active hyaluronidase compared to HYAL3. ³⁷ , ³⁸ These two transcripts were previously reported in the Johnson et al CUD GWAS using a different statistical method for TWAS (S‐PrediXcan). Additionally, a third gene TTC3 was identified using expression weights from the 3p21.3 putamen region (Z = 4.31, P = 1.63 × 10⁻⁵). Overexpression of TTC3 has been associated with negative effects on cognitive function and neuronal health. ³⁹ , ⁴⁰ As a result, this 3p21.3 region contains several plausible CUD risk genes, and future in silico and experimental follow‐up is warranted to disentangle true effects. In whole blood, ENO4 (Z = 4.8618, P = 3.26 × 10⁻⁶) and KIAA1598 (Z = 4.8618, P = 1.16 × 10⁻⁷), genes previously reported from European gene‐wise association analysis, were shown to be up‐regulated (Figure 1B). KIAA1598 is also known as SHTN1, with Shootin‐1 protein isoform switch from long isoform (Shtn1L) to short isoform (Shtn1S) known to play an important role in axonogenesis ⁴¹ , ⁴² (Table S4).

We expanded our analysis beyond mRNA expression by conducting a PWAS. In the brain, we again found that NAT6/HYAL1 was significantly up‐regulated in terms of protein abundance (Z = 5.019, P = 5.19 × 10⁻⁷) (Figure 1C). In the whole blood, we observed a novel significant association between predicted protein expression of NPTX1 and CUD (Z = −4.57, P = 4.88 × 10⁻⁶), which is critical for early human development ⁴³ (Figure 1D). Although NPTX1 down‐regulation was also detected in the TWAS, it did not reach statistical significance after correction for mRNA expression (Table S5).

The TWAS/PWAS analyses revealed HYAL1 as the most promising target for drug repurposing, as it can be modified by some existing nutraceuticals and investigational compounds (Table 1). However, the other TWAS/PWAS associations did not have current approved therapies that target them, although PDE4B (from the MAGMA analyses) is an existing drug target. If hyaluronidase dysregulation is indeed associated with CUD, as suggested by the TWAS/PWAS, it may represent a promising opportunity for therapeutic investigation. Hyaluronic acid, also known as hyaluronan, is a crucial component of the extracellular matrix and is known to play a role in the regulation of inflammatory processes and embryonic development. ⁴⁴ In the central nervous system (CNS), hyaluronic acid is a major component of perineuronal nets. These are mesh‐like structures that form around specific neuronal cell bodies and proximal dendrites and have a crucial role in synaptic stabilisation and plasticity. ⁴⁵ An increase in the expression of NAT6/HYAL1, which is associated with CUD, suggests that there may be abnormal breakdown of hyaluronic acid in individuals with this disorder, as hyaluronidases catalyse the degradation of hyaluronic acid. This abnormal breakdown could potentially increase the risk of developing CUD if the NAT6/HYAL1 gene has a direct causal effect. However, this interpretation is conflicted by the association we found between down‐regulation of another proximal hyaluronidase, HYAL3, and CUD. This requires further investigation, although there is evidence HYAL3 does not directly contribute to the metabolism of hyaluronic acid, ³⁷ but rather it is believed to do so by augmenting the activity of HYAL1. ³⁸ If hyaluronic acid catabolism were dysregulated in CUD, this finding may suggest the potential use of oral hyaluronic acid supplementation, which has been previously trialled for conditions like osteoarthritis. Additionally, the use of ascorbyl palmitate (L‐ascorbic acid 6‐hexadecanoate), an antioxidant that also inhibits hyaluronidase activity, may have therapeutic benefits as a neuroprotective agent. However, it is important to note that ascorbyl palmitate is not currently registered with the Food and Drug Administration (FDA) for therapeutic use (Table 1). ⁴⁶ , ⁴⁷ , ⁴⁸

To attempt to resolve the role of hyaluronic acid catabolism in CUD via effects at this locus, we applied probabilistic finemapping of TWAS test statistics on the region located at 3p21.3 that harbours both HYAL1 and HYAL3. Using a combined database of SNP weights from all GTEx tissues, as well some brain and blood datasets, HYAL1 was prioritised as the most likely causal gene for CUD, with a PIP of being in the 90% credible set of 89.3%, compared to 1.51% for HYAL3. However, probabilistic finemapping applied to the PsychENCODE SNP weight set, wherein HYAL1 did not have a model of genetically predicted expression available, prioritised HYAL3 instead (PIP = 78.3%). Given this disparity, further work is still required to reconcile how hyaluronidase activity may be involved in the aetiology of CUD.

We tested the genetic correlation between CUD and a panel of blood‐based biomarkers from the UKBB to gain further insight into drug repurposing opportunities by exploring the interplay between circulating biochemical factors and the pathophysiology of CUD. After Bonferroni correction, we found that CUD was genetically correlated with 12 of the 50 blood‐based biomarkers we tested. This included alanine aminotransferase (r _g  = 0.185, SE = 0.034, P = 7.25 × 10⁻⁷), CRP (r _g  = 0.206, SE = 0.55, P = 2.0 × 10⁻⁴), eosinophil count (r _g  = 0.122, SE = 0.032, P = 9.75 × 10⁻⁵), gamma glutamyltransferase (r _g  = 0.193, SE = 0.043, P = 7.94 × 10⁻⁶), lymphocyte count (r _g  = 0.178, SE = 0.033, P = 6.04 × 10⁻⁸), triglycerides (r _g  = 0.146, SE = 0.039, P = 2.0 × 10⁻⁴) and white blood cell (WBC) count (r _g  = 0.155, SE = 0.032, P = 1.89 × 10⁻⁶). It is worth nothing that inflammation of the CNS has long been linked to psychiatric disorders, including schizophrenia. ⁵⁷ Although we did not find direct evidence of a causal effect, the shared biology between CUD and the biomarkers we tested using LDSC is still informative for future treatment opportunities (Tables S6 and S7).

The authors confirm that the data supporting the findings of this study are publicly available and accessible online and can be found within the article and its supplementary material. The CUD summary statistics can be found here: 10.6084/m9.figshare.14842692.