Introduction
Substance use is a leading cause of global mortality (Ezzati, Lopez, Rodgers, Vander Hoorn, & Murray, Reference Ezzati, Lopez, Rodgers, Vander Hoorn and Murray2002) and shows a close association with psychiatric disorders. Individuals who have a prior lifetime history of psychiatric disorders (Conway, Swendsen, Husky, He, & Merikangas, Reference Conway, Swendsen, Husky, He and Merikangas2016) as well as those with a current diagnosis show higher rates of substance use across a variety of substances, including tobacco (Lasser et al., Reference Lasser, Boyd, Woolhandler, Himmelstein, McCormick and Bor2000; Weinberger et al., Reference Weinberger, Kashan, Shpigel, Esan, Taha, Lee and Goodwin2017), alcohol (Weitzman, Reference Weitzman2004), and cannabis (Blanco et al., Reference Blanco, Hasin, Wall, Flórez-Salamanca, Hoertel, Wang and Olfson2016). For example, relative to general population, individuals with severe psychotic disorders showed an increased risk of smoking [odds ratio (OR) = 4.6], heavy alcohol use (OR = 4.0), and heavy cannabis use (OR = 3.5) (Hartz et al., Reference Hartz, Pato, Medeiros, Cavazos-Rehg, Sobell, Knowles and Pato2014). Substance use has been found to predict psychiatric disorders prospectively, including depressive disorders (Brook, Cohen, & Brook, Reference Brook, Cohen and Brook1998), anxiety disorders (Johnson et al., Reference Johnson, Cohen, Pine, Klein, Kasen and Brook2000), and psychotic disorders (van Os et al., Reference van Os, Bak, Hanssen, Bijl, de Graaf and Verdoux2002). At the same time, early psychiatric problems have also been found to predict later substance use (King, Iacono, & McGue, Reference King, Iacono and McGue2004; Miettunen et al., Reference Miettunen, Murray, Jones, Mäki, Ebeling, Taanila and Moilanen2014). The elevated use of substances is associated with poorer outcomes such as medication non-adherence (Margolese, Malchy, Negrete, Tempier, & Gill, Reference Margolese, Malchy, Negrete, Tempier and Gill2004), more psychiatric hospitalizations (Dalton, Cate-Carter, Mundo, Parikh, & Kennedy, Reference Dalton, Cate-Carter, Mundo, Parikh and Kennedy2003), and increased suicide rates (Kask et al., Reference Kask, Ekselius, Brandt, Kollia, Ekbom and Papadopoulos2016), which is partly responsible for the excess mortality due to smoking and alcohol-related diseases (Brown, Inskip, & Barraclough, Reference Brown, Inskip and Barraclough2000; Hjorthøj et al., Reference Hjorthøj, Østergaard, Benros, Toftdahl, Erlangsen, Andersen and Nordentoft2015). The high comorbidity and accompanying mortality naturally leads to questions about shared etiology among these behaviors of high public health impact (Barkus & Murray, Reference Barkus and Murray2010; Gregg, Barrowclough, & Haddock, Reference Gregg, Barrowclough and Haddock2007).
Numerous studies have demonstrated shared genetic vulnerability to substance use and psychiatric disorder, especially those conceptualized as externalizing disorders such as ADHD and conduct disorder (Hicks, Krueger, Iacono, McGue, & Patrick, Reference Hicks, Krueger, Iacono, McGue and Patrick2004; Kendler et al., Reference Kendler, Aggen, Knudsen, Røysamb, Neale and Reichborn-Kjennerud2011; Rosenström et al., Reference Rosenström, Gjerde, Krueger, Aggen, Czajkowski, Gillespie and Ystrom2019; Young et al., Reference Young, Friedman, Miyake, Willcutt, Corley, Haberstick and Hewitt2009). Young et al. (Reference Young, Friedman, Miyake, Willcutt, Corley, Haberstick and Hewitt2009) reported a moderately heritable factor loading on substance use (alcohol, tobacco, marijuana, and other illicit drugs), conduct disorder, ADHD, and novelty-seeking traits at the start and the end of adolescence, using structural equation ACE models on the family data of 293 twin pairs. Similarly, Hicks et al. (Reference Hicks, Krueger, Iacono, McGue and Patrick2004) reported a highly heritable, general vulnerability factor that accounts for the correlations among conduct disorder, adult antisocial behavior, and alcohol and drug dependence in a sample of 542 twin families. Of note, most studies reporting such findings are based on known familial relationships (e.g. twins), rather than direct measurement of genetic markers. Genome-wide association studies (GWAS) have now identified hundreds of risk loci for substance use and psychiatric disorders, providing new opportunities to evaluate the pattern of relationships among these phenomena. Moreover, GWAS allows one to test not only the existence of shared genetic risk, but to interrogate individual genomic locations and related genes, with the potential to inform biology associated with the shared liability. The largest substance use GWAS to date identified many sets of genes and relevant tissues contributing to individual differences in alcohol and tobacco use such as nicotinic, dopaminergic, glutamatergic neurotransmission and developmental biology (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019). Gene-sets relevant to cannabis use also included neurogenesis and dopaminergic neurotransmission (Pasman et al., Reference Pasman, Verweij, Gerring, Stringer, Sanchez-Roige, Treur and Vink2018). Investigating genetic overlap can further elaborate knowledge on how these biological systems relate to psychiatric disorders or more specifically to certain types of substance use measures.
Pleiotropy – that variation within a given locus is associated with variation in two or more phenotypes – is abundant for complex traits (Paaby & Rockman, Reference Paaby and Rockman2013). Evidence for pleiotropy comes from family studies and population studies of genome-wide genetic covariances among phenotypes, as well as associations between multiple phenotypes and variation within individual genes or genomic loci. These avenues of research are informative about whether pleiotropy exists, but do not readily distinguish different types of pleiotropy. Horizontal pleiotropy arises when genetic variation affects two or more phenotypes independently [or, alternatively, independently affect some intermediate process(es), which then affects the two phenotypes]. Vertical pleiotropy arises when genetic variation affects one phenotype, the experience/expression of which then causes a second phenotype; this is the type of pleiotropy Mendelian randomization (MR) mainly attempts to detect. Genetic covariances can arise under either form of pleiotropy (Solovieff, Cotsapas, Lee, Purcell, & Smoller, Reference Solovieff, Cotsapas, Lee, Purcell and Smoller2013). Distinguishing horizontal from vertical pleiotropy would have implications for many existing causal theories of psychopathology and substance use, some of which posit that psychiatric disorders cause substance use (e.g. self-medication; Khantzian, Reference Khantzian1997), or that substance use is a causal risk factor of psychiatric disorders (Weiser & Noy, Reference Weiser and Noy2005). In contrast, the existence of horizontal pleiotropy would support conceptualizations of psychopathology as a product of more general liabilities to the manifestation of psychopathology and substance use or dependence (Kotov et al., Reference Kotov, Krueger, Watson, Achenbach, Althoff, Bagby and Zimmerman2017; Krueger & Markon, Reference Krueger and Markon2006). Such causal hypotheses can in theory be tested with MR through the use of genetic instruments (Davey Smith & Ebrahim, Reference Davey Smith and Ebrahim2003). Different MR methods provide complementary ways to test causal associations, and here we take advantage of two recently developed methods (O'Connor & Price, Reference O'Connor and Price2018; Pickrell et al., Reference Pickrell, Berisa, Liu, Ségurel, Tung and Hinds2016).
The aim of this paper is to inform the nature of shared genetic influences on substance use and psychiatric disorders by (1) examining the structure of genetic correlations using genome-wide methods and directly measured genetic variation (i.e. GWAS), (2) characterizing individual loci that are associated with two or more phenotypes, and (3) testing potential causal relations – horizontal or vertical pleiotropy – between substance use and psychiatric disorders. Although several GWAS meta-analyses have been conducted on the genetic overlap between various individual pairs of psychiatric disorders (Lee et al., Reference Lee, Anttila, Won, Feng, Rosenthal, Zhu and Smoller2019; van Hulzen et al., Reference van Hulzen, Scholz, Franke, Ripke, Klein, McQuillin and Reif2017), the overlap between multiple forms of substance use and a variety of psychiatric disorders has not been evaluated.
Materials and methods
Quality control
Publicly available European ancestry GWAS summary statistics, the largest to date, were collected for six substance use and five psychiatric disorder phenotypes. Any well-powered meta-analysis GWAS of substance use and psychiatric disorders covering externalizing, internalizing, and psychotic disorders were considered for inclusion. Among these, those reporting at least five significant loci in European ancestry and whose summary statistics were publicly available before mid-2019 were included. Additionally, education and height were included as controls to serve as a reference for interpreting the results. Quality controls, including variants filtering and LD clumping, were applied to ensure that problematic variants were excluded and different sets of summary statistics were as comparable as possible. Study characteristics after QC and details of QC procedure are reported in Table 1 and the Supplementary Note, respectively. Note that LDSC intercept values for the smoking and alcohol use traits are <1.0, consistent with the fact that these summary statistics were originally generated using genomic control (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019).
Table 1. Study information
a Number of GWAS variants are ascertained by LD clumping (r 2 > 0.1) genetic variants with p value <5 × 10−8 using Priority Pruner.
b Sample size is based on full summary statistics released in public except GSCAN phenotypes (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019) which included 23andMe sample.
c Smoking and drinking GWAS data from Liu et al. (Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019) have LDSC intercept <1 due to genomic control.
Factor analysis
We used genomicSEM (Grotzinger et al., Reference Grotzinger, Rhemtulla, de Vlaming, Ritchie, Mallard, Hill and Tucker-Drob2019) to perform factor analysis on the multivariate covariance matrix. The matrix was first constructed with SNPs filtered to those polymorphic in HapMap3 Europeans with minor allele frequency >0.01. We fit a hypothesized three-factor model as well as exploratory factor analysis (EFA) with 1–3 factors. Fit for all models were compared using confirmatory factor analysis (CFA). EFA was performed with promax rotation and maximum likelihood estimation using factanal() in the R package stats v.3.6.1. CFA was conducted with Diagonally Weighted Least Square estimation and implemented in genomicSEM, which provides standard fit information including the χ2 fit statistic, Akaike Information Criterion (Akaike, Reference Akaike1974; Vrieze, Reference Vrieze2012), comparative fit index (CFI), and standardized root mean square residual (SRMR). While no doubt the confirmatory models based on the EFA were overfitted, the CFA allowed us to evaluate fit indices for best-case (indeed, over-fitted) scenarios.
Bivariate locus-wise association
To identify pleiotropic loci, we used gwas-pw (Pickrell et al., Reference Pickrell, Berisa, Liu, Ségurel, Tung and Hinds2016; Ruderfer et al., Reference Ruderfer, Ripke, McQuillin, Boocock, Stahl, Pavlides and Kendler2018), a bivariate association method used to identify loci associated with pairs of phenotypes (here, 78 total pairs). The method is a Bayesian hierarchical modeling approach developed to compare the posterior probability of four competing models for a given approximately independent LD block. Models 1 and 2 assume the block harbors a causal variant that is associated with either only the first or only the second phenotype in a pair, respectively. Model 3 assumes the block harbors a variant that is associated with both phenotypes simultaneously, and model 4 assumes the block contains two distinct variants, each associated with only one of the two phenotypes. As in the original publication of the method, genomic regions were defined a priori by splitting chromosomes into approximately independent LD blocks (mean block size: 1.5 M base-pairs) (Berisa & Pickrell, Reference Berisa and Pickrell2016). The method empirically estimates priors for each regional model by using effect size information of variants in a given LD block and applies sample overlap correction (here we used effect size correlations between variants with p value >0.1 in both studies; online Supplementary Table S1) to generate posterior probabilities. A region was deemed pleiotropic if the posterior probability of model 3 was >0.9. SNPs showing the highest posterior probability within the shared locus were chosen as lead SNPs for that region. This method has been shown to detect shared loci at a false discovery rate of 10% (Pickrell et al., Reference Pickrell, Berisa, Liu, Ségurel, Tung and Hinds2016).
Functional annotations and tissue/gene-set enrichment analysis
First, we tested whether pleiotropic loci harbor more or less deleterious variants compared to non-pleiotropic loci using Combined Annotation Dependent Depletion (CADD) scores, which predict deleteriousness of variants based on protein function and structure (Kircher et al., Reference Kircher, Witten, Jain, O'Roak, Cooper and Shendure2014). Second, we tested whether variants from pleiotropic loci are more or less likely to harbor variants of different functional classes (e.g. non-synonymous, intronic, intergenic, etc.). Lead SNPs and their LD partners (r 2 > 0.4, ld-window 500 kb) were annotated with CADD and SEQMINER using refGene (retrieved 18 September 2019; Ye et al., Reference Ye, Krebs, Choukrallah, Keime, Plewniak, Davidson and Tora2010). Next, we performed tissue and gene-set enrichment analysis to test whether genes mapped to pleiotropic loci were enriched in certain tissue or gene-sets. A given region was categorized depending on whether it is associated with both substance use and psychiatric or associated with only one or the other. Tissue- and gene-set enrichment analyses were conducted by extracting nearest genes of the lead SNPs in these regions using a hypergeometric association test implemented in FUMA (Watanabe, Taskesen, van Bochoven, & Posthuma, Reference Watanabe, Taskesen, van Bochoven and Posthuma2017). Gene-sets surviving Bonferroni correction for tissue specificity p < 0.5/54, and FDR of 0.05 for gene-set analysis were considered significantly enriched.
Testing for vertical pleiotropy
To test causal associations, we applied bidirectional MR (Pickrell et al., Reference Pickrell, Berisa, Liu, Ségurel, Tung and Hinds2016) and the Latent Causal Variable model (LCV) (O'Connor & Price, Reference O'Connor and Price2018), two recently developed methods applicable to GWAS summary statistics. They are designed to be relatively robust to horizontal pleiotropy which is prevalent for complex traits (Verbanck, Chen, Neale, & Do, Reference Verbanck, Chen, Neale and Do2018). Both methods are based on simple intuition: if phenotype A causally influences phenotype B (e.g. smoking causing lung cancer), then any variant associated with A (smoking) will have correlated effects on B (lung cancer). However, variants associated with B (lung cancer) will not necessarily have correlated effects on A (smoking), as there will be many genetic causes of B (lung cancer) that are independent of A (e.g. asbestos exposure). In bidirectional MR, two sets of correlations are calculated on genetic association effect sizes, first using genome-wide significant (GWS) variants of the first phenotype and then repeating the same procedure using GWS variants for the second phenotype. These correlations are used to evaluate evidence for four models: in models 1 and 2, phenotype A or B causes B or A; in model 3, there are no causal relationships; in model 4, two phenotypes are very closely related (e.g. two alternate measurements of a single entity or one could be the major highly penetrant causal factor for the other). A relative likelihood for causal (models 1 and 2) v. non-causal (models 3 and 4) model (rl) <0.01 was interpreted as supporting the causality. In the LCV model, a latent variable mediates the genetic correlation between phenotypes A and B, having causal effects on both traits. A genetic causality proportion (GCP) is estimated (ranging from 0 to 1), which quantifies the degree of genetic causality, using mixed fourth moments of marginal effect sizes of all available SNPs. A high GCP indicates the genetic component of phenotype A is at least partially causal for phenotype B, while low GCP values suggest a lack of causal effects. Phenotype pairs with GCP >0.6 and p values surviving Bonferroni correction (here, p < 0.00064) were identified as potential phenotypes in causal relations (O'Connor & Price, Reference O'Connor and Price2018). Finally, the LCV and bidirectional MR results were compared to two-sample MR results obtained using MRbase v.0.5.3 (Hemani et al., Reference Hemani, Zheng, Elsworth, Wade, Haberland, Baird and Haycock2018) in R 3.6.1 on substance use–psychiatric phenotype pairs, for integration with the broad MR literature (see Supplementary Note for a detailed procedure for each method).
Results
Structure of genetic correlations
The first five eigenvalues of the genetic correlation matrix were 2.9, 1.1, 0.72, 0.13, and 0.07 (online Supplementary Fig. S1). We fit EFAs with 1–4 factors. Only the two-factor model was free of Heywood cases (loadings ⩾1.0), rendering model interpretation difficult. All EFA results are reported in online Supplementary Table S2. The two-factor model consisted of an externalizing factor loading on all three smoking phenotypes and ADHD, and an internalizing-psychosis factor loading on depression, anorexia, schizophrenia, and bipolar disorder (χ2 = 838.92, df = 26, p = 4.31 × 10−160, AIC = 876.92, CFI = 0.80, SRMR = 0.11). This procedure of fitting exploratory models and then similar confirmatory models in the same dataset will lead to overfitting, so these values should be considered an absolute best-case scenario. The four-factor model, roughly having an internalizing, psychosis, externalizing, and substance use factor, achieved the best fit. However, to converge, we had to drop drinks per week and constrain multiple uniquenesses to have >0 variance. Overall, models with 3–4 factors, including our hypothesized three-factor model, resulted in negative variances for several variables. As a sensitivity analysis, we conducted CFA using summary statistics for smoking and drinking to which no genomic controls had been applied, and this did not change the fit. CFA results are presented in online Supplementary Table S3 and Fig. S2.
Shared loci
Complete results on pleiotropic loci are available in online Supplementary Tables S4–8. Of all 30 possible substance use–psychiatric phenotype pairs, ever smoker and ADHD had the highest number of shared loci (N = 40), followed by age of smoking initiation and ADHD (N = 22). Within psychiatric disorder pairs, bipolar disorder and schizophrenia showed the highest number of shared loci (N = 108). Within substance use pairs, age of smoking initiation and ever smoker showed the highest number (N = 58). In contrast, height shared few loci with substance use and psychiatric disorders despite its large number of GWAS signals, but years of education showed substantial overlap with both substance use and psychopathology (Fig. 1). Of note, the number of pleiotropic regions reflects both degree of genetic overlap and statistical power of the original GWAS, rendering between-pair comparisons difficult. The lead SNPs in shared loci mostly agreed with the expected direction based on the genetic correlation (i.e. the more substance use, the higher risk of psychiatric disorder) (Fig. 1).
Fig. 1. The bottom triangle of panel A shows genetic correlations estimated by bivariate LDSC (online Supplementary Table S14). Non-significant correlations after Bonferroni correction are in gray. The upper triangle shows the number of loci shared between corresponding phenotype pairs. The numbers in the diagonal are the total number of genome-wide significant hits from original, univariate GWAS (ascertained by LD clumping). The background colors of the upper triangle indicate following pairs of domains: blue – between substance uses, purple – between substance use and psychiatric disorders, green – between psychiatric disorders, brown – between control phenotypes and all other phenotypes. Panel B shows the sign concordance of effect for lead SNPs in shared loci for a given phenotype pair. More reddish color indicates higher concordance. In each box, the left-side number indicates the number of lead SNPs having concordant directions of effects while the right-side number indicates the total number of lead SNPs in pleiotropic loci for a given pair. Note that smaller value of Age of Smoking Initiation indicates earlier age of smoking initiation and higher values of Smoking Cessation, Ever Smoker, and Lifetime Cannabis Use indicate current smoker, ever smoker, and ever cannabis user, respectively. Higher values of psychiatric disorders indicate the presence of the disorder.
A total of 341 loci showed evidence for association with two or more substance use or psychiatric disorder phenotypes. This number includes 230 loci reported in past GWAS of respective phenotypes, and 111 additional loci not reported in those studies (online Supplementary Table S8). Of these 341 loci, 152 loci (45%) were associated with both substance use and psychiatric phenotypes; 114 (33%) were associated with multiple substance use phenotypes but no psychiatric disorder, and 75 (22%) were associated with multiple psychiatric disorders but no substance use phenotypes .
Fig. 2. Panel A shows results for up-regulated Differentially Expressed Gene (DEG) set tests while panel B shows results from the down-regulated DEG set test. X-axis and Y-axis each represents tissue types and log10 p value. The bar higher than the dashed line indicates significantly enriched tissue type after Bonferroni correction. The white bar represents genes from loci pleiotropic for both substance use and psychiatric disorders while gray and black bars represent genes from loci pleiotropic within psychiatric disorders or substance use.
When 152 shared loci were classified depending on type of substance, 10 loci were associated with the use of all three types of substances [i.e. smoking, cannabis, and alcohol (Table 2)]. A gene encoding the D2 subtype of the dopamine receptor, DRD2, on chromosome 11 was located in these multi-substance shared loci along with cell-adhesion protein coding genes including NCAM1 and CADM2. Five loci were associated with alcohol and psychiatric disorder, including Utr3 region of FUT2, a gene involved in bifidobacterial diversity in the intestine and plasma levels of B12 vitamins (Mitchell, Conus, & Kaput, Reference Mitchell, Conus and Kaput2014). Five loci were associated specifically with cannabis use and psychiatric disorder, including intron region of WSCD2, a gene previously associated with extraversion and risky behaviors (Linnér et al., Reference Linnér, Biroli, Kong, Meddens, Wedow, Fontana and Nivard2019; Lo et al., Reference Lo, Hinds, Tung, Franz, Fan, Wang and Chen2017). Perhaps notably, nicotinic receptor genes (e.g. CHRNA4 and CHRNA5) did not belong to the 152 shared loci and instead were associated only with smoking (CHRNA4) and smoking and drinking phenotypes (CHRNA5).
Table 2. Pleiotropic regions associated with all three types of substance use and psychopathology
a Lead SNP in the locus, their associated phenotype pairs, and the direction of effects can be found in online Supplementary Tables S4 and 5.
b Two adjacent regions (3:84367614–85582078, 3:85582231–3:87409543) were combined to form this row.
Genes harboring lead SNPs of the 152 shared loci were overexpressed in all 13 brain tissues in GTEx v8 database and downregulated in six tissues including substantia nigra, kidney cortex, and liver (Aguet et al., Reference Aguet, Barbeira, Bonazzola, Brown, Castel, Jo and Lappalainen2019). They are enriched in 90 Gene Ontology (GO) biological processes, 10 molecular function, and 34 cellular components gene-sets. The top 3 most strongly associated biological processes were ‘neurogenesis’, ‘cell projection organization’, ‘regulation of nervous system development’. Enriched molecular functions include protein dimerization, transcription factor binding, and cell adhesion molecule binding, and the enriched cellular components encompass various parts of neuron including glutamatergic synapse. The 75 psychiatric-specific shared loci were overexpressed in four brain tissues (anterior cingulate cortex, frontal cortex, amygdala, and hippocampus) but not downregulated in any tissues. They were enriched in 11 biological processes, five molecular function, and 20 cellular components gene-sets. Most of the enriched gene-sets concerned neuronal signaling-related activity (e.g. voltage-gated ion channel activity) and cell structures (e.g. synaptic membrane). Finally, the 114 substance use-specific shared loci were overexpressed in the cerebellar hemisphere, cortex, and artery tibial, and downregulated in seven tissues including kidney cortex, liver, and pancreas. They are enriched in one molecular function, ‘clathrin binding’ and 12 cellular components most of which are located in neurons, including ‘dopaminergic synapse’ (Fig. 2; see online Supplementary Tables S9 and 10 for full tissue and gene-set enrichment results).
The median CADD rank scores were not significantly different across lead SNPs from substance-psychiatric pleiotropic loci, either psychiatric or substance use-specific pleiotropic loci, and non-pleiotropic loci (p = 0.34; online Supplementary Fig. S3). Most of the lead SNPs fell in intergenic or intronic regions (>90%) and their distribution across functional categories did not differ across different groups of pleiotropic and non-pleiotropic loci (p = 0.26; online Supplementary Fig. S4; online Supplementary Table S11).
Causally informative analyses
Of 78 pairs submitted to the analysis, none showed evidence of vertical pleiotropy from bidirectional MR and LCV. As a sensitivity analysis, bidirectional MR was conducted with the GWAS variants originally reported by each study and LCV was repeated with summary statistics for smoking and drinking without genomic control, both yielding the same pattern of results. Eight causal associations were significant in more than half of the four two-sample MR tests performed (Table 3). Full results of bidirectional MR and LCV, and two-sample MR are in online Supplementary Tables S12 and 13, respectively.
Table 3. Comparison of LCV, bidirectional MR and two-sample MR results
a The order of phenotypes in the first two rows represents the direction of causality (i.e. top to bottom).
b Age_Smoke = age of smoking initiation, Cannabis = lifetime cannabis use, MDD = major depressive disorder.
c Rho1, Rho2 = effect size correlations for genome-wide significant variants ascertained with phenotype 1 (Rho1) and phenotype 2 (Rho2).
d rl = relative likelihood (rl < 0.01 is considered significant causal association).
e Units: log odds ratio and standard deviations for binary and continuous phenotypes, respectively.
f Astericks attached to phenotype pairs represent significant causal associations after Bonferroni correction.
Discussion
We coordinated and analyzed GWAS results of six substance use and five psychiatric phenotypes to investigate the genetic correlation structure, multivariate association, and causal links among these phenotypes. The two-factor structure performed better than the single-factor model, but did not achieve a good model fit (e.g. CFI = 0.8, SRMR = 0.11, AIC = 876.92). Modification indices suggested that internalizing factor might not be well captured by covariance between depression and anorexia and residual covariance exists between schizophrenia and bipolar disorder as well as among some substance use phenotypes. In line with this diagnosis, the four-factor model achieved a better fit (e.g. CFI = 0.87, SRMR = 0.08, AIC = 709.62), which additionally fit separate substance use and psychosis factor. This model also aligns with current conceptualizations of the meta-structure of psychopathology (i.e. correlated externalizing, internalizing, and psychosis factors) (Kotov et al., Reference Kotov, Krueger, Watson, Achenbach, Althoff, Bagby and Zimmerman2017). However, caution is required in interpreting the four-factor model since the solution could be ill-specified (Heywood cases) and over-fitted (same data were used to both construct and test the model). The joint factor structure of substance use and psychiatric disorders should be further tested in independent samples, especially as large-scale GWAS on related phenotypes continue to be published.
We identified a total of 341 loci that showed at least one bivariate pleiotropic association with substance use or psychiatric phenotypes, which included about half of the loci reported in previous GWAS (~50% of the 458 univariate GWAS associations), confirming the presence of extensive pleiotropy. Genes nearest to the loci shared by both substance use and psychiatric disorders simultaneously (152 loci) were over-expressed in broad regions of the brain and enriched in neurodevelopmental pathways, suggesting that general neurodevelopmental processes (e.g. neurogenesis and neuron differentiation) may underlie the risk of both psychiatric disorders and substance use. For example, DCC, a gene involved in axonal growth and white matter projections (Jamuar et al., Reference Jamuar, Schmitz-Abe, D'Gama, Drottar, Chan, Peeva and Yu2017), was mapped to the most pleiotropic locus in a recent fixed-effect meta-analysis of eight psychiatric disorders (Lee et al., Reference Lee, Anttila, Won, Feng, Rosenthal, Zhu and Smoller2019). The same gene was mapped to highly pleiotropic locus in the current analysis, associated with the age of smoking initiation, ever smoker, depression, schizophrenia, and education. Biological adhesion and cellular morphogenesis/organization pathways were prominent among the neurodevelopmental pathways enriched in the 152 shared loci. Cell adhesion molecules have been reported to regulate synapse number, maturation, and plasticity (Sytnyk, Leshchyns'ka, & Schachner, Reference Sytnyk, Leshchyns'ka and Schachner2017), serving essential functions in neural development. On the other hand, genes nearest to the pleiotropic loci specific to psychiatric disorders were mostly enriched in neuronal signaling pathways such as voltage-gated ion channel activity while those specific to substance use were enriched in clathrin binding-related functions and cell structures. The potential role of clathrin binding is less clear, but it plays a key role in receptor endocytosis and vesicle recycling at th synapse (Jung & Haucke, Reference Jung and Haucke2007; Kaksonen & Roux, Reference Kaksonen and Roux2018) and has been implicated in drug-evoked neural plasticity in animal studies of addiction, including amphetamine (Brebner et al., Reference Brebner, Wong, Liu, Liu, Campsall, Gray and Wang2005), heroin (Van den Oever et al., Reference Van den Oever, Goriounova, Wan Li, Van der Schors, Binnekade, Schoffelmeer and De Vries2008), and morphine (Morón et al., Reference Morón, Abul-Husn, Rozenfeld, Dolios, Wang and Devi2007). For example, inhibiting clathrin-dependent postsynaptic AMPA receptor endocytosis in ventral mPFC and nucleus accumbens reduced drug-seeking behaviors in rats (Brebner et al., Reference Brebner, Wong, Liu, Liu, Campsall, Gray and Wang2005; Van den Oever et al., Reference Van den Oever, Goriounova, Wan Li, Van der Schors, Binnekade, Schoffelmeer and De Vries2008). Current results provide converging evidence as to the role of clathrin binding as a risk factor relatively specific to drug use and addiction.
Our bidirectional MR and LCV analysis detected no causally associated pairs, consistent with the notion that genetic correlations among these traits arise from horizontal pleiotropy. MR Egger, a classic two-sample method more robust to confounding by horizontal pleiotropy, also detected no pairs in vertical pleiotropy. In contrast, the other two-sample MR analyses more susceptible to such confounding (Verbanck et al., Reference Verbanck, Chen, Neale and Do2018) reported multiple pairs of causal associations, similar to past findings: genetic liability to ADHD associated with increased risk of smoking initiation (Fluharty, Sallis, & Munafò, Reference Fluharty, Sallis and Munafò2018) and vice versa (Treur et al., Reference Treur, Demontis, Smith, Sallis, Richardson, Wiers and Munafò2019), and genetic liability to smoking associated with increased risk of schizophrenia and depression and vice versa (Wootton et al., Reference Wootton, Richmond, Stuijfzand, Lawn, Sallis, Taylor and Munafòin press; Yao et al., Reference Yao, Xu, Cai, Liu, Ma, Li and Liin press); but see also some null MR findings (Gage et al., Reference Gage, Jones, Taylor, Burgess, Zammit and Munafò2017; Hodgson et al., Reference Hodgson, Coleman, Hagenaars, Purves, Glanville, Choi and Lewis2020; Taylor et al., Reference Taylor, Fluharty, Bjørngaard, Gabrielsen, Skorpen, Marioni and Munafò2014). Heterogeneity statistics from classic MR were significant for five out of eight associations (p < 0.05), indicating further the potential presence of horizontal pleiotropy for these five pairs (Bowden, Hemani, & Davey Smith, Reference Bowden, Hemani and Davey Smith2018). Three pairs, i.e. ADHD lowering the age of smoking initiation, ever smoker increasing risk for ADHD, and depression, showed non-significant heterogeneity statistics thus may be causally related. However, their genetic correlations were rather high (i.e. rg = −0.6, 0.57, 0.35 for each pair in order), a condition which can produce excess false positives in classic MR (O'Connor & Price, Reference O'Connor and Price2018).
On the other hand, LCV and bidirectional MR showed relatively well-calibrated type 1 and 2 errors in such a scenario (O'Connor & Price, Reference O'Connor and Price2018). Taken together, vertical pleiotropy, at least a single causal direction, may not explain the well-known high rate of co-occurrence among psychopathology and substance use. Related causal hypotheses such as self-medication or substance use as a causal risk factor for psychopathology are not well supported by current results. For example, GCP for major depression and ever smoker (0.04, s.e.: 0.06) and that for cannabis use and schizophrenia (−0.15, s.e.: 0.20) were close to zero and significantly different from 1.
Our findings should be interpreted in light of several limitations. First, the number of manifest variables may be still insufficient to fit CFA models with more than three factors. We also stress again that CFA fit could be overestimated to some unknown degree since the models were built on EFA results in the same sample. Current results should be tested in replication sample with more manifest variables especially those for internalizing. From a more technical standpoint, some of the constituent GSCAN cohorts used linear mixed models (LMM) and this may have exerted a subtle influence on CFA fit due to model misspecification between linear regression-based LDSC model and LMM (Yang, Zaitlen, Goddard, Visscher, & Price, Reference Yang, Zaitlen, Goddard, Visscher and Price2014). The degree and remedy for this potential bias should be further studied. For pleiotropic loci analysis, the size of LD blocks used in this method is relatively large (mean 1.5 Mb) yielding limited resolution, and the posterior probability threshold to detect pleiotropy is almost certainly less stringent than the typical family-wise error rate of 5% in GWAS. We used the more liberal threshold (corresponding to an FDR of ~10%) since the focus of the current study was to assess the scope of pleiotropy in known associated loci, rather than declaring high confidence novel loci. In the simulation performed in the original article, gwas-pw overall gave a slight overestimation for model 3 (pleiotropy). In the presence of high sample overlap, it additionally gave modest overestimation for model 4 (two separate effects in a single locus) over model 1 or 2 (only one effect for either trait in a single locus), requiring caution in interpretation. It is also not possible to distinguish a single variant influencing both traits from two causal variants in strong LD affecting each trait separately. Regarding causal inference, although the results were mainly interpreted in light of shared genetic effects, cyclical feedback loops between phenotypes may exist to some degree, possibly for all pairs of phenotypes. This type of causality can evade LCV and bidirectional MR but can still be detected by two-sample MR methods to some extent. Last, the GWAS were all performed on individuals of European ancestry, which limits the generalizability of current findings. Increasing availability of trans-ancestry GWAS can improve the generalizability of this kind of analyses in the future.
Despite these limitations, the current study provides comprehensive analyses of the genetic overlap between substance use measures and psychiatric disorders based on GWAS results. It offers novel insights into the structure, biology, and causal nature of this overlap and points toward the future directions for methodological and etiology research. Refining our understanding of the biological underpinnings of the externalizing spectrum (e.g. salient developmental processes and periods) will be helpful to tackle adverse outcomes associated with it. Continuing to study multiple layers of biology such as those underlying addiction to various substance types (e.g. clathrin-mediated endocytocis of AMPA receptors) and that pertaining to specific substance–psychiatric disorder association can add rich information to the etiology of the two phenomena. Finally, the current study applied novel causal inference methods developed to address pervasive horizontal pleiotropy in complex traits and adds to the fast-growing causal inference literature.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S003329172000272X.
Acknowledgements
The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing resources that contributed to the research results reported within this paper. We would like to thank the research participants and employees of 23andMe for making this work possible. Members of the 23andMe research team are Michelle Agee, Babak Alipanahi, Adam Auton, Robert K. Bell, Katarzyna Bryc, Sarah L. Elson, Pierre Fontanillas, Nicholas A. Furlotte, David A. Hinds, Karen E. Huber, Aaron Kleinman, Nadia K. Litterman, Matthew H. McIntyre, Joanna L. Mountain, Elizabeth S. Noblin, Carrie A.M. Northover, Steven J. Pitts, J. Fah Sathirapongsasuti, Olga V. Sazonova, Janie F. Shelton, Suyash Shringarpure, Chao Tian, Joyce Y. Tung, Vladimir Vacic, and Catherine H. Wilson. We would also like to thank Hoang Nguyen for his advice on factor analysis.
Financial support
This study was supported by funding from US National Institutes of Health awards R01DA037904 and R01DA044283 to S.V.; R01HG008983 and R21DA040177 to D.J.L.
Conflict of interest
The authors associated with the 23andMe Research Team are employees of 23andMe Inc. All other authors reported no biomedical financial interests or potential conflicts of interest.
Ethical standards
Ethical review and approval were provided by the University of Minnesota institutional review board.
