Introduction
Prevalence of smoking in schizophrenia (SCZ) patients is considerably larger than in general population (Sagud, Mihaljevic Peles, & Pivac, Reference Sagud, Mihaljevic Peles and Pivac2019). Therefore, if not taken into account, this may have some consequences in the characterization of the genetic risk for SCZ. For instance, two associated regions including genes coding nicotinic acetylcholine receptor subunits, specifically, the CHRNA2 gene and the CHRNA5–CHRNA3–CHRNB4 cluster, were associated with SCZ and several smoking behaviors in large genome-wide association studies (GWAS) (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019; Ripke et al., Reference Ripke, Neale, Corvin, Walters, Farh, Holmans and O'Donovan2014). CHRNA2 gene was also associated with cannabis use disorder (CUD) (Demontis et al., Reference Demontis, Rajagopal, Thorgeirsson, Als, Grove, Leppälä and Børglum2019). Although our previous work revealed that the main association signals for SCZ and smoking at the CHRNA5–CHRNA3–CHRNB4 cluster were independent, it demonstrated that the signal at the CHRNA2 gene colocalized (Al-Soufi & Costas, Reference Al-Soufi and Costas2021). This may be due to a causal effect of substance use on SCZ (i.e. vertical pleiotropy) or to the fact that the same genetic factors may confer risk to SCZ and substance use independently (i.e. horizontal pleiotropy) (Solovieff, Cotsapas, Lee, Purcell, & Smoller, Reference Solovieff, Cotsapas, Lee, Purcell and Smoller2013). Bidirectional causality is also possible, for instance, if the substance causes risk to SCZ (Gage & Munafò, Reference Gage and Munafò2015), and if SCZ patients use substance to reduce symptoms (Goff, Henderson, & Amico, Reference Goff, Henderson and Amico1992; Kumari & Postma, Reference Kumari and Postma2005).
To know if the association of these and other genes with SCZ was driven by substance use, the analysis should be done by sample stratification based on substance use (Gage & Munafò, Reference Gage and Munafò2015). Previously, we found that a polygenic risk score on SmkInit was associated with SCZ in a case–control study, but the association was lost after inclusion of smoking status as a covariate (Al-Soufi et al., Reference Al-Soufi, Martorell, Moltó, González-Peñas, García-Portilla, Arrojo and Costas2022).
Mendelian randomization (MR) studies analyze putative causal association of exposure with outcome based on the use of SNPs associated with an exposure as instrumental variables, i.e., proxies for the exposure (Davies, Holmes, & Davey Smith, Reference Davies, Holmes and Davey Smith2018; Treur, Munafò, Logtenberg, Wiers, & Verweij, Reference Treur, Munafò, Logtenberg, Wiers and Verweij2021). SNPs have to accomplish three key assumptions to be valid instrumental variables for MR studies: (1) they are robustly associated with exposure (the relevance assumption), (2) there are no confounders of the association between SNPs and outcome (the independence assumption), and (3) SNPs affect the outcome only through their effect on the exposure (the exclusion restriction assumption). Using this methodology, an effect of smoking on SCZ has been described (Barkhuizen, Dudbridge, & Ronald, Reference Barkhuizen, Dudbridge and Ronald2021; Wootton et al., Reference Wootton, Richmond, Stuijfzand, Lawn, Sallis, Taylor and Munafò2020; Yuan, Yao, & Larsson, Reference Yuan, Yao and Larsson2020), in agreement with epidemiological data (Gurillo, Jauhar, Murray, & MacCabe, Reference Gurillo, Jauhar, Murray and MacCabe2015; Kendler, Lönn, Sundquist, & Sundquist, Reference Kendler, Lönn, Sundquist and Sundquist2015; Koyanagi, Stickley, & Haro, Reference Koyanagi, Stickley and Haro2016), although the results were not confirmed by all sensitivity analyses. In addition, one relevant method, multi-trait-based conditional and joint analysis (mtCOJO), allows for conditioning a GWAS summary statistics on the GWAS summary statistics for the putative causal trait (Byrne et al., Reference Byrne, Zhu, Qi, Skene, Bryois, Pardinas and Wray2021; Zhu et al., Reference Zhu, Zheng, Zhang, Wu, Trzaskowski, Maier and Yang2018). This method has been used to analyze the effect of genetic susceptibility to major depression disorder (MDD) on suicide GWAS (Mullins et al., Reference Mullins, Kang, Campos, Coleman, Edwards, Galfalvy and Willour2022), the effect of SCZ on breast cancer (Byrne et al., Reference Byrne, Ferreira, Xue, Lindström, Jiang, Yang and Chenevix-Trench2019), or the identification of disorder-specific risk factors after adjustment for general psychopathology (Byrne et al., Reference Byrne, Zhu, Qi, Skene, Bryois, Pardinas and Wray2021). Furthermore, change in genetic correlation with other traits after conditioning may reveal the impact of the adjustment factor on the shared susceptibility (Hill et al., Reference Hill, Davies, Ritchie, Skene, Bryois, Bell and Deary2019; Mullins et al., Reference Mullins, Kang, Campos, Coleman, Edwards, Galfalvy and Willour2022).
The GWAS and Sequencing Consortium of Alcohol and Nicotine use (GSCAN) found a significant genetic correlation between smoking behavior and both externalizing traits, such as number of drinks per week (DrnkWk), lifetime cannabis use, attention deficit-hyperactivity disorder (ADHD), risk tolerance, or age when having the first child (AFB); as well as internalizing traits (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019), such as MDD or the related personality trait of neuroticism (Griffith et al., Reference Griffith, Zinbarg, Craske, Mineka, Rose, Waters and Sutton2010). This agrees with the existence of two alternative pathways to substance use escalation and dependence: an externalizing pathway, related to novelty-seeking, impulsivity and difficulties in self-regulation, and an internalizing pathway, characterized by substance use to cope with negative emotional states (Dir & Hulvershorn, Reference Dir and Hulvershorn2019; Green et al., Reference Green, Conway, Silveira, Kasza, Cohn, Cummings and Compton2018; Heitzeg, Hardee, & Beltz, Reference Heitzeg, Hardee and Beltz2018). There was also a significant genetic correlation between smoking behavior and SCZ. How this shared genetic susceptibility with SCZ is related to internalizing and externalizing phenotypes is currently unknown.
In this work, under MR assumptions, we applied mtCOJO to condition the largest GWAS summary statistics of SCZ from the Psychiatric Genomics Consortium Wave 3 (PGC-SCZ3) (Trubetskoy et al., Reference Trub, Pardiñas, Qi, Panagiotaropoulou, Awasthi, Bigdeli and O'Donovan2022) on smoking with the following main goals, (i) to identify loci that lost significance, indicative of their association with SCZ being at least partially mediated by smoking; (ii) to identify new significant loci whose effect on SCZ was hidden due to confounding by smoking, (iii) to characterize anatomical regions, temporal states and gene-sets differentiating SCZ genetic susceptibility from that shared with smoking, and (iv) to identify genetic correlations among SCZ and substance use traits, psychiatric disorders, behavioral traits and neuropsychological traits affected by smoking. Additionally, colocalization analysis was performed to identify specific loci that confirm the general findings.
Methods
An overview of the study design is shown in online Supplementary Fig. S1.
GWAS summary statistics
Summary statistics for PGC-SCZ3, restricted to European ancestry cohorts, was used as the main data (Trubetskoy et al., Reference Trub, Pardiñas, Qi, Panagiotaropoulou, Awasthi, Bigdeli and O'Donovan2022). Four different smoking behavior GWAS from the GSCAN were also used: smoking initiation, i.e., ever v. never regular smokers (SmkInit), smoking cessation, i.e., current v. former smokers (SmkCes), number of cigarettes smoked per day, either as a current smoker or as a former smoker (CigDay) and age started smoking regularly (AgeSmk) (Liu et al., Reference Liu, Jiang, Wedow, Li, Brazel, Chen and Vrieze2019) (Table 1).
Table 1. Main GWAS summary statistics and GSMR analysis results using each smoking trait as the exposure trait and SCZ as outcome trait
a Some of the original GWAS included 23andMe samples, not used in the present analysis.
b Number of independent significant SNPs used at GSMR analysis after exclusion of pleiotropic SNPs.
c Number of pleiotropic SNPs detected by HEIDI-outlier test.
d As estimated in Kahn et al. (Reference Kahn, Sommer, Murray, Meyer-Lindenberg, Weinberger, Cannon and Insel2015).
e As estimated in Bobes, Arango, Garcia-Garcia, and Rejas (Reference Bobes, Arango, Garcia-Garcia and Rejas2010).
f As estimated in Creamer et al. (Reference Creamer, Wang, Babb, Cullen, Day, Willis and Neff2019).
GWAS summary statistics for different traits of interest comprising substance use traits, psychiatric disorders, behavioral traits, and neuropsychological traits were used for genetic correlation analysis (online Supplementary Table S1). These GWAS were selected according to their large sample size, public availability, and expected or detected genetic correlation with SCZ.
Conditional analysis
In order to select the smoking trait with the highest causal effect on SCZ to be used in conditional analysis by mtCOJO, generalized summary data-based Mendelian randomization (GSMR) was used to estimate the causal effect ($\hat{b}_{xy}$
) of each smoking trait (x) from GSCAN, as exposure trait, on the outcome trait SCZ (y) (Zhu et al., Reference Zhu, Zheng, Zhang, Wu, Trzaskowski, Maier and Yang2018), as implemented in GCTA (Yang, Lee, Goddard, & Visscher, Reference Yang, Lee, Goddard and Visscher2011). GSMR uses SNPs at r 2 < 0.05 and p < 5 × 10−8 for exposure trait as genetic instruments to estimate $\hat{b}_{xy}$
. HEIDI-outlier test was used to remove SNPs with pleiotropic effects (Zhu et al., Reference Zhu, Zheng, Zhang, Wu, Trzaskowski, Maier and Yang2018). This test is based on detection of heterogeneity in $\hat{b}_{xy}$
estimated at the different SNPs used as instrumental variables under the idea that the expected value of $\hat{b}_{xy}$
estimated at any of the SNPs is identical in absence of pleiotropy. All parameters for GSMR analysis were set as default, except that a minimum number of SNPs, set as 10 by default, was removed. In order to get interpretable estimates of $\hat{b}_{xy}$
for binary traits, GSMR estimates were transformed to the liability scale, as described in Byrne et al. (Reference Byrne, Zhu, Qi, Skene, Bryois, Pardinas and Wray2021) (online Supplementary Methods).
mtCOJO, as implemented in GCTA (Yang et al., Reference Yang, Lee, Goddard and Visscher2011), was used to condition the effect of SNPs at the PGC-SCZ3 on the selected smoking GWAS of GSCAN (Byrne et al., Reference Byrne, Zhu, Qi, Skene, Bryois, Pardinas and Wray2021; Zhu et al., Reference Zhu, Zheng, Zhang, Wu, Trzaskowski, Maier and Yang2018). Specifically, mtCOJO estimates the effect of each SNP on SCZ conditioning on the effect of smoking on SCZ (i.e. $\hat{b}_{xy}$
estimated by GSMR) as: $\hat{b}_{zy}$
|$\hat{b}_{xy}$
= $\hat{b}_{zx}$
– $\hat{b}_{zx} \hat{b}_{xy}$
, where $\hat{b}_{zx}$
and $\hat{b}_{zy}$
are the estimated effects of each SNP z on the exposure x and the outcome y, respectively. mtCOJO analysis is robust to sample overlap thanks to the incorporation in the model of the sampling covariance between SNP effects. mtCOJO also performs Linkage Disequilibrium Score regression (LDSC) to obtain genetic correlation, SNP-based heritability, and intercept of the bivariate LDSC. Linkage Disequilibrium (LD) Scores files based on The 1000 Genomes phase 3 European samples (T1000G_Eur) were used for this analysis. Sample prevalence and population prevalence of each binary trait were supplied to mtCOJO in order to calculate SNP-based heritability on the liability scale (Table 1). T1000G_Eur were also used as LD reference sample to get independent SNPs for GSMR and to remove SNPs with a large difference of allele frequency among GWAS summary-statistics and this LD reference sample. mtCOJO estimates conditional effect size, conditional standard error, and conditional p values for each shared SNP between GWAS. Simulation studies indicated that mtCOJO is robust to collider bias and the estimates of $\hat{b}_{xy}$
by GSMR are almost identical to those from MR Egger-regression, an alternative method considered to be free of confounding by pleiotropy (Zhu et al., Reference Zhu, Zheng, Zhang, Wu, Trzaskowski, Maier and Yang2018).
Significant loci definition
The SNP2GENE module of FUMA was used to get significant independent loci at r 2 < 0.1, or more than 250 kb apart, and p < 5 × 10−8 and to annotate, for each locus, all significant SNPs at p < 0.05 in LD with the most significant SNP, i.e., the lead SNP (Watanabe, Taskesen, Van Bochoven, & Posthuma, Reference Watanabe, Taskesen, Van Bochoven and Posthuma2017). The T1000G_Eur was used as the reference panel for LD. SNPs from reference panel not included in GWAS were also considered for annotation. Major histocompatibility complex (MHC) was not excluded to define significant loci.
Overlapping significant loci between original SCZ GWAS and conditional analysis were defined as those loci for which the lead SNP at one analysis is in LD at r 2 > 0.1 with the lead SNP at the other analysis. Lost loci were those significant loci at original SCZ GWAS that were no longer significant at conditional analysis (conditional p > 5 × 10−8). New loci were those significant loci at conditional analysis but not at original SCZ GWAS (conditional p < 5 × 10−8).
Loci were defined as concordant or discordant for two traits according to the effect of the most significant shared SNP, i.e., concordant if this SNP is associated with both traits in the same direction and discordant if this SNP is associated with both traits in opposite directions.
Heritability and genetic correlation analysis
LDSC was used to estimate SNP-based heritabilities and genetic correlations from GWAS summary statistics (Bulik-Sullivan et al., Reference Bulik-Sullivan, Finucane, Anttila, Gusev, Day, Loh and Neale2015). LD Scores based on T1000G_Eur were used and GWAS summary statistics were filtered to HapMap3 SNPs. Summary statistic filters were set as default. SNP-based heritability was transformed to the liability scale using sample prevalence and population prevalence (Table 1). A significantly different genetic correlation between the original GWAS and the conditional analysis with the trait of interest was assessed using a two-sided test [2 × pnorm(− abs(abs(r gi − r gj)/sqrt(s.e.i2 + s.e.j2)))], as in Hill et al. (Reference Hill, Davies, Ritchie, Skene, Bryois, Bell and Deary2019).
Colocalization analysis
Colocalization analysis was performed following Pickrell et al. (Reference Pickrell, Berisa, Liu, Ségurel, Tung and Hinds2016) in order to identify if GWAS association signals for two traits at the same locus are due to the same causal variant (see online Supplementary Methods for further details).
Gene mapping
Gene mapping of significant loci, restricted to protein coding genes, was done with the SNP2GENE module of FUMA in several steps (Watanabe et al., Reference Watanabe, Taskesen, Van Bochoven and Posthuma2017). First, genes were selected by position, considering genes up to 10 kb apart from the SNPs annotated at each significant locus. Second, genes for which annotated SNPs were cis-eQTLs at a FDR < 0.05 in brain tissues were selected, using the following databases: Brainseq, PsychENCODE, xQTLServer, CommonMind Consortium, BRAINEAC and GTEx v8 for brain. The rest of parameters were set as default. Finally, a prioritization step was done based on functional mapping. Specifically, we selected the most significant annotated SNP with a putative functional consequence on a gene or genes, considering exonic non-synonymous or frameshift substitution SNPs, SNPs affecting splicing, SNPs within an UTR or a cis-eQTLs, and prioritized the corresponding gene or genes. If the most significant SNP is a cis-eQTL for several genes, we prioritized the gene for which that SNP is the most significant eQTL. The protein coding gene nearest to the lead SNP was also indicated.
Enrichment analysis
Enrichment analysis was performed with MAGMA using the SNP2GENE module in FUMA to test for enrichment for association at our conditional analysis and the original SCZ GWAS in differentially expressed genes in tissues (using data from GTEx v8), brain developmental stages and ages (data from BrainSpan), and precompiled Gene Ontology (GO) gene-sets (c5.bp, c5.cc and c5.mf) from MsigDB (Watanabe et al., Reference Watanabe, Taskesen, Van Bochoven and Posthuma2017). MHC region was excluded. We assigned to each gene those SNPs within 35 kb upstream and 10 kb downstream, as previously done (O'dushlaine et al., Reference O'dushlaine, Rossin, Lee, Duncan, Parikshak, Newhouse and Breen2015; Pardiñas et al., Reference Pardiñas, Holmans, Pocklington, Escott-Price, Ripke, Carrera and Walters2018; Trubetskoy et al., Reference Trub, Pardiñas, Qi, Panagiotaropoulou, Awasthi, Bigdeli and O'Donovan2022).
Results
Conditional analysis
The causal effect of each of the smoking traits on SCZ was assessed using GSMR. SmkInit and AgeSmk were the smoking traits with a higher effect on SCZ in the liability scale (Table 1). However, estimation using AgeSmk was based only on 7 instrumental SNPs, lacking reliability. Based on the largest sample size, the largest number of instrumental SNPs used in GSMR analysis and the most significant effect on SCZ, SmkInit was chosen as the exposure trait to condition SCZ GWAS.
SNP-based heritability of SCZ after conditioning on SmkInit (SCZcond) was 0.1901 (s.e. = 0.0060), similar to the original GWAS (0.1998, s.e. = 0.0069). Genetic correlation between original SCZ GWAS and this conditional analysis was high (r g = 0.9726, s.e. = 0.0012) but significantly different from 1 (p = 2.14 × 10−115). Genetic correlation between SCZcond and SmkInit significantly changed in comparison with original SCZ GWAS (p = 8.55 × 10−15) and became negative (r g = −0.0764, s.e. = 0.0206, p = 2.00 × 10−4).
Conditioning SCZ GWAS on SmkInit led to a reduction of the number of independent significant loci from 174 at original SCZ GWAS to 150 at SCZcond (Fig. 1). One significant locus at the original SCZ GWAS, chr8:3506797–5059268, was not included in the conditional analysis because there were no SNPs in LD with the lead SNP at this region in the SmkInit GWAS.
Fig. 1. Miami plot comparing original SCZ GWAS and SCZcond. The x axis shows chromosomal position, and the y axis shows the significance of association. Red line indicates the genome-wide significance level (p = 5 × 10−8). Lost loci at conditional analysis are highlighted in blue and new loci in orange. For each lost and new locus, prioritized gene and nearest protein coding gene to the lead SNP are shown. Loci for which colocalization analysis supports a common association signal between SCZ and SmkInit are indicated by *.
Overlapping loci between SCZ and SCZcond
From the 174 significant loci at original SCZ GWAS, 131 loci remained significant at SCZcond (online Supplementary Table S2). Effect size of these overlapping loci did not significantly change after conditioning (0.3791% decrease, paired t test: p = 0.5318) indicating that they were primarily loci with an effect on SCZ independent from SmkInit.
Lost loci in SCZcond
Forty-two of the loci associated with SCZ at original GWAS lost association after conditioning on SmkInit, suggesting that their association with SCZ may be partially driven by smoking, under MR assumptions (Fig. 1 and online Supplementary Table S3). Unexpectedly, although almost all loci were concordant between SCZ and SmkInit, three of them were discordant. However, they had a higher effect on SCZcond suggesting that they lost significance due to inappropriate fit of the LD reference panel.
Colocalization analysis supported a common association signal between SCZ and SmkInit for eight of these lost loci (Fig. 1). Six of them were genome-wide significant for SmkInit.
New loci in SCZcond
Conditioning on SmkInit led to the identification of 19 new significant loci for SCZ that, under MR assumptions, may not have been detected in the original SCZ GWAS due to the confounding effect of smoking (Fig. 1 and online Supplementary Table S4). Colocalization between SmkInit and SCZcond was supported for two of these new loci, encompassing MRPL19/GCFC2 and LRRTM4 (Fig. 1). Twelve of the 19 new loci were reported as significant loci in the extended PGC3 GWAS based on multiple ancestries and additional replication samples, highlighting the validity of the analysis to identify new associated loci (online Supplementary Table S4).
Enrichment analysis before and after conditioning on SmkInit
Tissue enrichment analysis revealed a similar pattern for the original SCZ GWAS and SCZcond (Fig. 2a and online Supplementary Fig. 2a). Regarding developmental stages, prenatal stages were more associated with SCZcond while later stages were more enriched in the original SCZ GWAS, suggesting that smoking genetic susceptibility has an effect on temporal enrichment pattern for SCZ (Fig. 2b). Early mid-prenatal stage was the significant stage common to both GWAS results. Late mid-prenatal was only significantly enriched in SCZcond and late infancy was only significant for SCZ before conditioning. Using more specific brain ages did not lead to any clear difference in enrichment between the original SCZ GWAS and SCZcond, suggesting that enrichment was based on genes expressed at more broad temporal periods (online Supplementary Fig. S2b). Brain ages enrichment was more significant for SCZ before conditioning.
Fig. 2. MAGMA enrichment analysis with original SCZ GWAS and SCZcond. The solid vertical line represents the significance at the Bonferroni's correction for the number of gene expression profiles tested. a: significantly enriched tissues, b: brain developmental stages enrichment.
Many GO terms related to synapsis were significant. In general, GO terms were more associated with the original SCZ GWAS than with SCZcond. Some of them reached significance only in SCZ. All but one GO term enriched for SCZcond, the cellular component ‘membrane protein complex’, were also enriched at the original SCZ GWAS (online Supplementary Fig. S3).
Changes in genetic correlation of SCZ with other traits of interest after conditioning on SmkInit
Genetic correlation between SCZ, before and after conditioning, and several traits was assessed in order to analyze the impact of smoking on the shared genetic susceptibility between SCZ and other phenotypes, under MR assumptions (Fig. 3 and online Supplementary Table S5). Genetic correlation of SCZ with some smoking phenotypes significantly changed after conditioning on SmkInit (Fig. 3a). Genetic correlation with lifetime smoking index (LifetimeSmk) and SmkCes was no longer significant; while negative genetic correlation with AgeSmk became positive. Regarding other substance use traits, correlation with lifetime cannabis use was significantly reduced after conditioning, while genetic correlation with DrnkWk was completely lost. For other traits, there was a non-significant change in the expected direction.
Fig. 3. Genetic correlation (r g) between different traits (a: substance use traits, b: substance use disorders and psychiatric disorders, c: other traits) and SCZ before and after conditioning on SmkInit. Traits are described in online Supplementary Table S1. Error bars represent ± 1 s.e.. Genetic correlations with original SCZ GWAS and with SCZcond are indicated by a circle and a triangle, respectively. Genetic correlations that are significant at a FDR < 0.05 are colored in black. A significantly different genetic correlation between original SCZ GWAS and SCZcond with the trait of interest using a two-sided test is indicated by: *: p < 0.05, **: FDR < 0.05. LifetimeSmk, lifetime smoking index; SmkCes, smoking cessation; Cigday, number of cigarettes smoked per day; AgeSmk, age started smoking regularly; DrnkWk, number of drinks per week; AUDIT, Alcohol Use Disorder Identification Test Total Score; TUD, tobacco use disorder; CUD, cannabis use disorder; Alcohol dep, alcohol dependence; Opioid dep, opioid dependence; ADHD, attention deficit-hyperactivity disorder; MDD, major depressive disorder; BIP, bipolar disorder; ASD, autism spectrum disorder; AN, anorexia nervosa; OCD, obsessive-compulsive disorder; AD, anxiety disorders; Risk PC1, first PC of the four risky behaviors in the UK Biobank; AFB, age when having the first child; AFS, age at first sexual intercourse; SWB, subjective well being; EA, years of educational attainment.
Regarding substance use disorders (SUD), changes in the expected direction were observed in all cases, reaching significance for tobacco use disorder (TUD) and CUD (Fig. 3b). Among non-substance-related psychiatric disorders, conditioning SCZ on SmkInit led to a significant change in correlation just for ADHD.
Some behavioral traits were also differentially correlated with SCZ and SCZcond (Fig. 3c). Genetic correlation of SCZ was significantly reduced after conditioning for Townsend Social Deprivation Index, number of sexual partners and the first PC of the four risky behaviors in the UK Biobank (Risk PC1), including SmkInit (online Supplementary Table S1). Likewise, negative genetic correlation of SCZ with AFB and age at first sexual intercourse (AFS) became positive after conditioning. Finally, genetic correlation of SCZ with years of educational attainment (EA) became significant and positive after conditioning.
Colocalization analysis with other traits of interest
Colocalization analysis was done for the 14 traits with a significantly different genetic correlation with SCZ before and after conditioning to strengthen the evidence of shared genetic susceptibility. Among the 42 lost loci at SCZcond, 25 presented some evidence of association with other trait (p < 1 × 10−5 for a SNP in LD r 2 > 0.1 with the SCZ lead SNP), involving all the analyzed traits except TUD and social deprivation (online Supplementary Table S6). Colocalization analysis supported a common association signal between the original SCZ GWAS and the other trait for 14 of the 25 loci (Fig. 4). As expected, these loci had concordant effects between SCZ and those traits whose genetic correlation with SCZ was reduced after conditioning on SmkInit. Similarly, the shared lost loci were discordant between SCZ and those traits whose genetic correlation with SCZ increased after conditioning on SmkInit. The only exception was the locus that maps to CHRNA2/CLU, whose association signal for SCZ showed concordant effects with LifetimeSmk and CUD but was discordant for the number of sexual partners. A more complicated pattern was found for EA. As expected by the increase in the genetic correlation between SCZ and EA after conditioning, three of the four lost loci that colocalized between SCZ and EA were associated in opposite directions with both traits; while one of them was concordant.
Fig. 4. Colocalization for lost loci shared with a trait with a significantly different genetic correlation with SCZ before and after conditioning. Only those loci with evidence of colocalization are shown. The x axis shows the significance of the most significant SNP for SCZ which is also significant at p < 1 × 10−5 for the trait of interest. For each locus, prioritized gene and nearest protein coding gene to the SCZ lead SNP are shown. Those loci for which colocalization analysis also supports a common association signal between SCZ and SmkInit are indicated with an asterisk. Loci for which the shared SNP has concordant effects between SCZ and the trait of interest are indicated in bold type. Direction of GWAS effect estimates was reversed for AgeSmk, AFB and AFS to indicate association with a higher externalizing behavior. LifetimeSmk, lifetime smoking index; AgeSmk, age started smoking regularly; DrnkWk, number of drinks per week; CUD, cannabis use disorder; ADHD, attention deficit-hyperactivity disorder; Risk PC1, first PC of the four risky behaviors in the UK Biobank; AFB, age when having the first child; AFS, age at first sexual intercourse; EA, years of educational attainment.
Regarding the 19 new loci at SCZcond, there was some evidence of association for nine of them with six of the analyzed traits. Colocalization analysis did not support a common association signal between SCZcond and the corresponding trait for any of these loci, either because of a lack of power or because they were independently associated with each trait (online Supplementary Table S7).
Discussion
Here we assessed the effect of the shared genetic susceptibility between smoking behavior and SCZ on the largest SCZ GWAS using mtCOJO. Under MR assumptions, we identified new putative risk loci for SCZ whose association may be reduced by smoking acting as a confounding factor and loci whose association with SCZ at the original GWAS may be partially driven by smoking. Several specific examples of colocalization of the association signal at new and lost loci between SCZ and SmkInit provided additional evidence that the change in significance was due to the effect of these loci on susceptibility to smoking. Brain developmental stages enrichment analysis found that prenatal expression is more important in SCZ after adjustment by SmkInit, suggesting a difference in temporal expression of susceptibility. Conditioning SCZ on smoking resulted in a significant change of the genetic correlation with substance use and dependence traits, ADHD, and several behavioral traits.
In addition, our results are consistent with previous MR findings showing evidence for a causal effect of several smoking traits on SCZ (Wootton et al., Reference Wootton, Richmond, Stuijfzand, Lawn, Sallis, Taylor and Munafò2020) using a different method, GSMR. GSMR presents some advantages, such as the detection of pleiotropic SNPs as HEIDI outliers and the consideration of residual LD between instrumental variables. A GSMR analysis was previously performed to estimate the effect of SmkInit and CigDay on SCZ, using previous SCZ GWAS, obtaining similar results (Barkhuizen et al., Reference Barkhuizen, Dudbridge and Ronald2021; Byrne et al., Reference Byrne, Ferreira, Xue, Lindström, Jiang, Yang and Chenevix-Trench2019). Interestingly, we have found a significant, although small, negative genetic correlation between SmkInit and SCZ after conditioning on SmkInit. This result is similar to that of a recent study reporting that SmkInit was negatively genetically correlated with SCZ after accounting for the genetic effects of cannabis use, CUD and nicotine dependence by genomic structural equation modeling (Johnson et al., Reference Johnson, Hatoum, Deak, Polimanti, Murray, Edenberg and Agrawal2021).
Our genetic correlation analysis revealed a significant reduction in genetic correlation between externalizing phenotypes, such as substance use and SUD or ADHD, and SCZ after conditioning on SmkInit. By contrast, genetic correlation with SCZ after conditioning on SmkInit did not change for internalizing phenotypes, such as MDD, anxiety disorder or neuroticism. This is a clear indication that the shared genetic susceptibility between SCZ and SmkInit is mainly related to the externalizing pathway towards tobacco use. This agrees with the higher manifestation of impulsivity found in schizophrenic patients in comparison with healthy controls (Heerey, Robinson, McMahon, & Gold, Reference Heerey, Robinson, McMahon and Gold2007; Nolan, D'Angelo, & Hoptman, Reference Nolan, D'Angelo and Hoptman2011; Ouzir, Reference Ouzir2013). As in general population, high levels of impulsivity in schizophrenic patients may lead to substance use (Dervaux et al., Reference Dervaux, Baylé, Laqueille, Bourdel, Le Borgne, Olié and Krebs2001; Gut-Fayand et al., Reference Gut-Fayand, Dervaux, Olié, Lôo, Poirier and Krebs2001). Nevertheless, this association between SCZ and impulsivity was not replicated in other studies suggesting heterogeneity in SCZ or differences in the analyzed dimension of impulsivity (Fischer et al., Reference Fischer, McMahon, Kelly, Wehring, Meyer, Feldman and Gorelick2015; Reddy et al., Reference Reddy, Lee, Davis, Altshuler, Glahn, Miklowitz and Green2014).
There were a few traits whose genetic correlation with SCZ after conditioning on SmkInit did not fit the interpretation based on externalizing liability. On the one hand, automobile speed propensity and several measures of irritability may be considered externalizing traits, but genetic correlation did not change. This agrees with the poor fit of a genomic structural equation model testing the existence of a common latent factor among these traits and other externalizing traits (Karlsson Linnér et al., Reference Karlsson Linnér, Mallard, Barr, Sanchez-Roige, Madole, Driver and Dick2021). On the other hand, EA and social deprivation showed significant differences in genetic correlation with SCZ after conditioning. Both traits are genetically correlated with a putative latent externalizing factor based on seven different externalizing traits (Karlsson Linnér et al., Reference Karlsson Linnér, Mallard, Barr, Sanchez-Roige, Madole, Driver and Dick2021). The interpretation of the genetic basis of EA is very complex, and assortative mating on correlated traits, such as social deprivation, plays a role (Okbay et al., Reference Okbay, Wu, Wang, Jayashankar, Bennett, Nehzati and Young2022).
Using mtCOJO, we detected some loci that reached significance in the conditional analysis and other that lost significance. However, the p values did not dramatically change in any case. Although some changes in significance may be simply caused by random fluctuations of p value due to estimation errors, this may be alternatively explained by the existence of a direct effect of these SNPs on SCZ, in addition to their effect on SCZ driven by SmkInit. Some of these variants are associated with SCZ and SmkInit in the same direction, leading to a loss of significance when conditioning on SmkInit, and others in the opposite direction, resulting in a gain of significance after conditioning. This agrees with the recent evidence of extensive pleiotropy with both same and opposite direction among psychiatric disorders (Hindley et al., Reference Hindley, Frei, Shadrin, Cheng, O'Connell, Icick and Andreassen2022). Interestingly, 12 of the 19 new loci identified in the European GWAS after conditioning on SmkInit were also significant in the extended GWAS that used mixed ancestry and additional replication samples (Trubetskoy et al., Reference Trub, Pardiñas, Qi, Panagiotaropoulou, Awasthi, Bigdeli and O'Donovan2022), confirming the validity of our approach.
Our colocalization analysis confirmed the existence of several pleiotropic genes affecting externalizing traits and SCZ among those genes that lost significance after conditioning. Interestingly, some of these loci have been previously identified in multitrait analyses. Thus, eight loci encompassing RERE, ARTN/PTPRF, CHRNA2/CLU, RNF122/DUSP26, DDHD2/LETM2, SORCS3, BCL11B and DCC were significantly associated with a cross-disorder meta-analysis of the PGC, including eight psychiatric disorders without inclusion of SUD (Lee et al., Reference Lee, Anttila, Won, Feng, Rosenthal, Zhu and Smoller2019). Similarly, two loci encompassing NAT8/ALMS1 and PCDH7 were significantly associated with a latent externalizing factor based on seven traits (Karlsson Linnér et al., Reference Karlsson Linnér, Mallard, Barr, Sanchez-Roige, Madole, Driver and Dick2021). Finally, five loci encompassing ARTN/PTPRF, CUL3, DDHD2/LETM2, SORCS3, and METTL15 were previously highlighted as loci jointly associated with SCZ and risky behaviors (Hindley et al., Reference Hindley, Bahrami, Steen, O'Connell, Frei, Shadrin and Andreassen2021).
In the case of LSAMP, which colocalized with several externalizing traits in our analysis, animal models confirmed its role in social behavior, anxiety, hyperactivity and heightened response to environmental stressors (Catania, Pimenta, & Levitt, Reference Catania, Pimenta and Levitt2008; Philips et al., Reference Philips, Lilleväli, Heinla, Luuk, Hundahl, Kongi and Vasar2015).
One locus that deserves special attention is CHRNA2. At this gene, we confirmed the previously found common association signal between SCZ, smoking, and CUD (Al-Soufi & Costas, Reference Al-Soufi and Costas2021). We also found colocalization at this locus with the number of sexual partners, but in the opposite direction to that expected if externalizing behavior was involved. So, these findings support an implication of this gene in substance use that may be partially driving its association with SCZ regardless of externalizing susceptibility.
Two of the new risk loci for SCZ, encompassing CDH11 and TOM1L2/GID4, were significant at the previous SCZ GWAS from the PGC (Ripke et al., Reference Ripke, Neale, Corvin, Walters, Farh, Holmans and O'Donovan2014). This result may be indicative of a higher proportion of smokers in the latest SCZ GWAS, highlighting the existence of genetic heterogeneity in patients with SCZ (Alnæs et al., Reference Alnæs, Kaufmann, Van Der Meer, Córdova-Palomera, Rokicki, Moberget and Westlye2019; Jablensky, Reference Jablensky2006).
The main limitations of the study are those associated with the GWAS used in the analyses. Many of them, mainly those for behavioral traits, are based on ‘minimal phenotyping’ approach (Cai et al., Reference Cai, Revez, Adams, Andlauer, Breen, Byrne and Flint2020). In addition, GWAS capture not only direct effects of the genotype on the studied trait, but also indirect effects such as genotype-environment correlations, uncorrected population stratification, or assortative mating on their trait or correlated traits (Okbay et al., Reference Okbay, Wu, Wang, Jayashankar, Bennett, Nehzati and Young2022; Plomin & von Stumm, Reference Plomin and von Stumm2022), which may affect interpretation of results. Moreover, only those publicly available GWAS were used, at expenses of not always having the largest sample size, which limits the statistical power of the different analyses. Another limitation is related to colocalization analysis. Success to identify colocalization may be reduced in case of sample overlap or the existence of several causal variants at the locus. However, colocalization was performed to confirm general genomic findings. Thus, although we probably did not detect all colocalized loci, the impact for the study conclusions is minimal. Finally, we must take into account the limitations of the MR approach, which relies on several assumptions. Thus, although our GSMR analysis supports a causal role of several smoking behaviors on SCZ risk, we cannot rule out this effect being driven by cannabis smoking (Reed, Wootton, & Munafò, Reference Reed, Wootton and Munafò2022). This may be due to vertical pleiotropy, with both smoking and cannabis being on the same path to SCZ. In this case, the MR exclusion restriction assumption is not violated. But it may also be a case of horizontal pleiotropy with SNPs influencing both traits independently. Other pathways that violate this assumption may also exist, such as maternal smoking during pregnancy. In fact, there is a clear association between maternal smoking during pregnancy and SCZ in offspring (Hunter, Murray, Asher, & Leonardi-Bee, Reference Hunter, Murray, Asher and Leonardi-Bee2020; Niemelä et al., Reference Niemelä, Sourander, Surcel, Hinkka-Yli-Salomäki, McKeague, Cheslack-Postava and Brown2016). However, a large population- and family-based epidemiological study revealed that the more likely explanation for this association is familiar confounding rather than a true causal effect on risk (Quinn et al., Reference Quinn, Rickert, Weibull, Johansson, Lichtenstein, Almqvist and D'Onofrio2017). Due to the limitations of MR, triangulation of evidence using different approaches is needed to firmly establish a role of smoking on SCZ risk (Lawlor, Tilling, & Smith, Reference Lawlor, Tilling and Smith2016; Ohlsson & Kendler, Reference Ohlsson and Kendler2020). Strong triangulation evidence is accumulating in recent years that smoking negatively affects mental health (Firth, Wootton, Sawyer, & Taylor, Reference Firth, Wootton, Sawyer and Taylor2023).
In summary, we applied mtCOJO to condition SCZ GWAS on smoking behavior, leading to the identification of potential new risk loci for SCZ whose effects may have been masked by smoking, risk loci partially associated to SCZ through smoking, and evidence of shared genetic susceptibility between SCZ and smoking behavior specifically related to externalizing behavior. Taking into account the high prevalence of smoking in patients with psychiatric disorders, as well as some evidence for causality (Firth et al., Reference Firth, Wootton, Sawyer and Taylor2023; Treur et al., Reference Treur, Munafò, Logtenberg, Wiers and Verweij2021), it may be useful to extend our approach to other psychiatric disorders to increase our knowledge of the role of smoking, a modifiable risk factor, on mental health. Extension to other substances may also be useful, although, in general, there is currently a limitation due to the scarcity of SNPs robustly associated with the exposure.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291723000326
Acknowledgements
The authors wish to thank all the Consortia and research groups providing the summary statistics datasets used in this work.
Financial support
This work was supported by the Instituto de Salud Carlos III (grant number ISCIII/PI17/01424, cofounded by FEDER, and grant number RIAPAd RD21/0009/0011 to JC). LA-S was supported by a research fellowship from the Ministerio de Universidades (grant number FPU18/03243).
Conflict of interest
None.
