Cannabis use and the risk of primary open-angle glaucoma: a Mendelian randomization study

Study design

MR uses genetic variants as instrumental variables to assess causal associations between risk factors and diseases based on the random assignment of genetic variants in individuals at conception⁷. These genetic variants are usually single-nucleotide polymorphisms (SNPs) and since they are randomly allocated in individuals independently of other factors, MR studies can serve as naturally occurring randomized controlled trials⁷. Thus, MR association estimates are less prone to biases occurring from confounding and reverse causation than those derived from traditional observational studies. We conducted a two-sample, summary-based MR and utilized summary statistics from three genome-wide association studies (GWAS) of lifetime cannabis use⁸, cannabis use disorder⁹ and POAG¹⁰. Then by combining these estimates the causal association between cannabis use and cannabis use disorder with POAG was calculated. The recommendations by STROBE-MR¹¹ and “Guidelines for performing Mendelian randomization investigations” were followed⁷. The study protocol was not pre-registered.

Data sources

We retrieved summary data from the largest GWAS to date for lifetime cannabis use comprising 184,765 individuals of European descent, by the International Cannabis Consortium, 23andMe, and UK Biobank⁸. The exposure was defined as any self-reported use of cannabis during a person’s lifetime. GWAS analysis were adjusted for sex, age, ancestry, and genotype batch. Genotyping and imputation methods have been described elsewhere⁸. We also retrieved summary statistics for cannabis use disorder from a GWAS meta-analysis of 17,068 cases and 357,219 controls of European descent, derived from the Psychiatric Genomics Consortium Substance Use Disorders working group, Lundbeck Foundation Initiative for Integrative Psychiatric Research (iPSYCH), and deCODE (Supplementary Table S1)⁹. Cases from the Psychiatric Genomics Consortium had the diagnosis of cannabis abuse or dependence according to the Diagnostic and Statistical Manual of Mental Disorders (DSM)-IV or DSM-III-R, from clinician ratings or semi-structured interviews. IPSYCH cases met the criteria for a diagnosis of cannabis abuse (F12.1) or cannabis dependence (F12.2) based on the ICD-10 criteria, while cases from the deCODE sample were diagnosed with lifetime cannabis abuse or dependence according to DSM-IV or DSM-III-R, or with cannabis use disorder according to DSM-V. Genotyping, quality control and imputation methods have been described elsewhere⁹. SNP-POAG associations were taken from a GWAS meta-analysis of 16,677 POAG cases and 199,580 controls of European ancestry¹⁰ from 16 participating studies (Supplementary Table S1). POAG was defined according to ICD9/ICD10 criteria. GWAS adjusted for age, sex, and study-specific principal components¹⁰. Genotyping, quality control and imputation have been described in detail elsewhere¹⁰. There was a 12.7% overlap between the GWAS of lifetime cannabis use and POAG, which does not significantly affect our association estimates.

Selection of genetic variants as instrumental variables

We adopted two approaches in the selection of genetic variants as instrumental variables. In the primary analysis we selected only SNPs reaching genome-wide significance (P-value < 5*10^–8 for lifetime cannabis use and P-value < 5*10^–7 for cannabis use disorder) following clumping for linkage disequilibrium (LD) at r² < 0.001 across a 10mb window. In our secondary, more liberal approach^12,13, we selected SNPs independently associated with lifetime cannabis use and cannabis use disorder at a GWAS P-value < 5*10^–5 after accounting for LD at r² < 0.1, in order to increase pooled instrument strength and power of the analysis. In both approaches, we calculated the percentage of phenotypic variance that is explained by our exposures of interest. By summing the coefficients of determination (R²) derived from the associations of the selected SNPs with our exposures of interest, we were able to calculate the proportion of variability in our exposure phenotypes that is explained by the selected SNPs. Finally, we performed the MR-Steiger directionality test to identify the direction of causality between lifetime cannabis use and POAG and removed SNPs that were more strongly correlated with the outcome than the exposure¹⁴. We excluded SNPs with highly influential data points in the funnel plots and scatter plots of SNP–exposure and outcome associations. Five SNPs associated with lifetime cannabis use and eleven SNPs associated with cannabis use disorder were selected in the primary analysis, while 267 and 157 SNPs associated with lifetime cannabis and cannabis use disorder, respectively, were selected in the secondary analysis.

Statistical analysis

After data harmonization, where SNPs were filtered according to HapMap3¹⁵, excluded if they were strand-ambiguous and their effect sizes were aligned, we calculated Wald ratios by dividing the per-allele logarithm of odds ratio (logOR) for each selected SNP from the lifetime cannabis use and cannabis use disorder GWAS by the corresponding logOR from the same SNP in the POAG GWAS. Then, we estimated the effect of lifetime cannabis use and cannabis use disorder on the risk of POAG by pooling the Wald ratios with multiplicative random effects inverse-variance weighted (IVW) meta-analyses¹².

Univariable two-sample MR was performed using summary-level statistics from the largest available GWAS on lifetime cannabis use and POAG. The two-sample MR approach rests on 3 core assumptions: (1) the genetic instruments should be robustly associated with the exposure of interest (“relevance” assumption), (2) the genetic instruments are not associated with confounders of the exposure-outcome association (“exchangeability” assumption), and (3) the genetic instruments are associated with the outcome exclusively through their effect on the exposure of interest (“exclusion restriction” assumption)^16,17. The “relevance” assumption is satisfied by selecting SNPs, as instrumental variables, reaching the genome-wide significance (P-value < 5*10^–8). Moreover, in order to quantify instrument strength, we calculated the proportion of variance of the exposure explained by the genetic instruments, as well as the F-statistic of our instruments¹⁸. Although, the “relevance” and “exclusion restriction” assumption cannot be proven, we performed sensitivity analyses to assess any possible violations of these assumptions. These can occur through horizontal pleiotropy, where the genetic variants affect the outcome via biological pathways other than the exposure under investigation. Thus, in the primary analysis, we utilized PhenoScanner¹⁹ to investigate associations between the selected genetic instruments with traits that could potentially confound our analysis and in case that pleotropic pathways were discovered, multivariable MR was used to adjust for these effects²⁰. More specifically, one of our instrumental SNPs for lifetime cannabis use was associated with previously reported obesity-related phenotypes (Supplementary Table S3). Several studies have found an association between body mass index and POAG^21,22, so this SNP might have been associated with POAG through pathways other than our exposure of interest and, thus, we performed multivariable IVW adjusting for BMI. Additionally, the associations of each selected SNP and its proxies (r² > 0.8) with known risk factors for POAG were also checked. In the multivariable MR analyses the conditional F-statistic was used as a quantification of the strength of our genetic instruments²³. Moreover, we assessed the heterogeneity among the selected genetic variants in the primary analysis through the Cochran Q heterogeneity test and I_GX²¹⁷ in order to detect pleiotropy. MR Egger regression was performed in order to assess the presence of directional pleiotropy¹⁷, as well as pleiotropy-robust methods²⁴ (penalized weighted median, IVW radial regression and MR-Pleiotropy Residual Sum and Outlier (MR-PRESSO)). Since only five SNPs for lifetime cannabis use were selected in our primary analysis, the IVW radial regression and the MR-PRESSO were not performed²⁴. In order to assess whether the IVW estimate was driven by a single SNP, leave-one-out analysis was also conducted.

In our secondary analysis using a liberal threshold, we performed multiplicative random-effects IVW and pleiotropy-robust methods (penalized weighted median, IVW radial regression, MR-Pleiotropy Residual Sum and Outlier (MR-PRESSO))²⁴. The CAUSE MR analysis was additionally conducted as an additional method to improve statistical power and mitigate the risk of weak instrument bias^7,25.

All MR estimates for the associations between our exposures and POAG were multiplied by log_e2 (= 0.693), representing the change in log odds of POAG per doubling in the prevalence of our exposures²⁶. All analyses were performed with R version 4.2.1²⁷ using the MendelianRandomization, TwoSampleMR, MVMR, MR-PRESSO and cause packages.