Anti-Epileptic Drug Target Perturbation and Intracranial Aneurysm Risk: Mendelian Randomization and Colocalization Study

Background: In a genome-wide association study of intracranial aneurysms (IA), enrichment was found between genes associated with IA and genes encoding targets of effective anti-epileptic drugs. Our aim was to assess if this pleiotropy is driven by shared disease mechanisms that could potentially highlight a treatment strategy for IA. Methods: Using 2-sample inverse-variance weighted Mendelian randomization and genetic colocalization analyses we assessed: (1) if epilepsy liability in general affects IA risk, and (2) whether changes in gene- and protein-expression levels of anti-epileptic drug targets in blood and arterial tissue may causally affect IA risk. Results: We found no overall effect of epilepsy liability on IA. Expression of 21 genes and 13 proteins corresponding to anti-epileptic drug targets supported a causal effect (P<0.05) on IA risk. Of those genes and proteins, genetic variants affecting CNNM2 levels showed strong evidence for colocalization with IA risk (posterior probability>70%). Higher CNNM2 levels in arterial tissue were associated with increased IA risk (odds ratio, 3.02; [95% CI, 2.32–3.94]; P=3.39×10−16). CNNM2 expression was best proxied by rs11191580. The magnitude of the effect of this variant was greater than would be expected if systemic blood pressure was the sole IA-causing mechanism in this locus. Conclusions: CNNM2 is a driver of the pleiotropy between IA and anti-epileptic drug targets. Administration of the anti-epileptic drugs phenytoin, valproic acid, or carbamazepine may be expected to decrease CNNM2 levels and therefore subsequently decrease IA risk. CNNM2 is therefore an important target to investigate further for its role in the pathogenesis of IA.

The inverse variance weighted (IVW) MR method was selected as our main MR analysis method. Two linkage disequilibrium (LD)-independent variants are needed to perform IVW. IVW estimates are accurate if all genetic variants included conform to the MR assumptions of no heterogeneity. If only a single LD-independent variant was available, we used a Wald ratio test to quantify its MR effect but did not use it for further analysis since with a single genetic variant no test for pleiotropy can be performed, making it impossible to assess reliability of the estimates. MR-PRESSO was used in addition to exclude eQTL and pQTL where pleiotropic outliers influenced the IVW estimate. The MR-PRESSO method is an extension of the IVW method and is robust for horizontal pleiotropy amongst the genetic variants compared to IVW, by identifying and removing pleiotropic outliers from the analysis. 45,46 First, MR-PRESSO was used to identify pleiotropic outliers, which were then removed, and MR-PRESSO was subsequently used to calculate the MR effect using R package MRPRESSO (version 1.0). Empirical p-value was obtained by 1000 permutations. At least four LDindependent variants are needed to perform MR-PRESSO.
Instruments (i.e., uncorrelated variants) for MR should be associated with the exposure. To achieve this, a threshold for association p-value is typically applied. The downside of a stringent threshold is that it limits the number of variants available thereby potentially decreasing accuracy. This means that no single optimal p-value threshold exists. Therefore, we used two p-value thresholds for the association of a variant with the exposure: p < 5×10 -5 and p < 5×10 -8 . After this selection the genetic variants were clumped to select LD-independent instruments (LD r 2 > 0.01, or outside 10 megabase window to include) using the 1000 Genomes European reference panel. The analysis was performed using the R packages TwoSampleMR (version 0.5.5) and MRInstruments (version 0.3.2). 43

Colocalization analysis
If two traits show a genetic association signal on the same genomic region, colocalization analyses can estimate the probability that these signals are driven by the same causal variant. The colocalization approach assumes that 1) there is a maximum of one causal genetic instrument for either trait, 2) the causal probability of a genetic instrument is independent of the causal probability of other genetic instruments in the analysis, and 3) every causal genetic variant (genotyped or imputed) is included in the colocalization analysis. Following these assumptions, there are five hypotheses for each performed analysis: Hypothesis 0 (H0) states there are no causal genetic instruments present for either trait; (H1) states there is one causal genetic instrument present for the trait 1; (H2) states there is one causal genetic instrument present for trait 2; (H3) states both traits have a causal variant in the region, but these are distinct; (H4) states there is one shared causal genetic instrument present for both traits. 25 COLOC requires two prior probabilities. The first prior defines the probability tan any variant is causal for any trait. This was set to 1×10 -4 (i.e. one in 10,000 genome-wide genetic instruments are causal for either trait). The second prior defines the probability that any variant is causal for both traits and was set to 1 x 10 -6 . 25

Sensitivity analyses
The estimation of a causal effect by Mendelian randomization will be biased if the majority of the genetic instruments are invalid or if any of the genetic variants show pleiotropic effects. Therefore, the MR Egger regression, weighted median-based, weighted mode-based, simple mode-based and Mendelian Randomization Pleiotropy RESidual Sum and Outlier (MR-PRESSO) robust methods were performed in extension of the inverse variance weighted method. The claim of any causal effect of the exposure on intracranial aneurysms will be more credible when the results of the MR analyses are consistent. 46 The MR Egger regression and MR-PRESSO methods detect pleiotropic effects. MR Egger estimates the causal effect from the weighted regression slope, without forcing the intercept through zero. The results of the MR Egger method will be uniform with the effect size estimate of the IVW if the MR-Egger intercept shows a value of zero. This method requires independent pleiotropic effects of the genetic variants associated with the exposure, which allows for every genetic variant used in the MR analysis to have any form of pleiotropic effect. 47 The weakness of the MR Egger method is that the estimation is specifically affected by genetic variant outliers with large effects, and violation of the InSIDE assumption (instrument strength of direct effect, meaning effect of genetic variants should be independent of direct effects of genetic variants on the outcome). 48 The weighted and simple mode methods require that most genetic variants contribute to the estimation of the true causal effect instead of other occurring effects. So, both weighted-and modebased methods are more robust to a small number of pleotropic variants compared to the MR-Egger and IVW MR methods. 49,50 The Cochran's Q test was done to quantify the heterogeneity in the effects of genetic variants. If heterogeneity is large, this could indicate that Mendelian randomization assumptions are violated. Cochran's Q-test required the genetic variants to be independent. 51 The Steiger test of directionality is used to infer the causal direction of exposure and outcome. It is robust to measurement error, even when the effect of exposure on outcome is small. MR Steiger requires an assumed direction and then tests whether this direction is consistent with a causal effect. It indicated a causal relationship (in any direction) by testing the hypothesis that the difference in correlations between genetic variant g with exposure x and g with outcome y is greater than expected under the null hypothesis of identical correlations. The direction is inferred from the MR analysis and Steiger test. 52

Assessing the role of blood pressure in the CNNM2 locus
We tested whether the effect of increased in CNNM2 gene expression, genetically proxied by the most likely causal variant rs11191580, on IA, could be explained solely by increased blood pressure caused by rs11191580. If this is true, then the MR effect size estimate of blood pressure genetically proxied by only rs1119580 (the exposure) on IA (the outcome) is similar to the MR effect size estimate of blood pressure genetically proxied by all other genetic variants associated with blood pressure combined (the exposure) on IA (the outcome). Therefore, we performed a Wald ratio test of blood pressure (systolic and diastolic) on IA using only variant rs1119580 and an IVW test of genome-wide variants associated with blood pressure (P-value threshold = 5×10 -5 ) on IA excluding variants in the CNNM2 gene +/-1 megabases.
To further validate that the effect that genetically proxied differential CNNM2 levels have on IA is not solely mediated through blood pressure, we performed a conditional analysis. We conditioned the IA summary statistics on summary statistics for blood pressure, similar to using blood pressure as covariate in the IA genome-wide association study. Conditioning was done using multi-trait conditional and joint (mtCOJO) analysis 53 using its default setting using summary statistics for systolic and diastolic blood pressure measurements from the UK Biobank together (http://www.nealelab.is/uk-biobank/). European-ancestry genotypes from the 1000 Genomes project were used as reference panel.

Benchmarking the role of CNNM2 in blood pressure and IA
Since we limited our primary analyses to anti-epileptic drug targets and thus other genes with a stronger effect on IA that we did not include may be present, we benchmarked the importance of CNNM2 in IA and blood pressure to all other genes. We conducted MR with expression of each gene in tibial artery, and in whole blood as exposures, and blood pressure, and IA as outcomes. To allow a fair comparison of effect sizes, we standardized genome-wide association effect size of gene expression to mean=0, variance=1 for each gene separately. Then, we performed an inverse variance weighted MR analysis of expression of each gene separately on IA, and on systolic and diastolic blood pressure in both tibial artery and blood.

Supplementary Fig. S2. Gene-wise comparison of Mendelian randomization (MR) effects on blood pressure and on intracranial aneurysms in tibial artery.
Each point represents a gene of which genetically proxied expression has an MR effect on A) diastolic blood pressure or B) systolic blood pressure with MR P-value < 0.001. Error bars and shaded area: 95% confidence intervals. CNNM2 is highlighted. Lines are based on a linear regression of the per-gene MR effects of IA on blood pressure including all genes (red) and excluding CNNM2. Figure on next page.

Supplementary Fig. S3. Gene-wise comparison of Mendelian randomization (MR) effects on blood pressure and on intracranial aneurysms in blood.
Each point represents a gene of which genetically proxied expression has an MR effect on A) diastolic blood pressure or B) systolic blood pressure with P-value < 0.001. Error bars and shaded area: 95% confidence intervals. CNNM2 is highlighted. Lines are based on a linear regression of the per-gene MR effects of IA on blood pressure including all genes (red) and excluding CNNM2. Figure on next page.

Section
Checklist item    Figure S1 shows scatterplots of regional association statistics of intracranial aneurysms, CNNM2 expression, and IRF1 expression, to support our colocalization analyses.