Characterization of migraine risk loci

We mapped the 123 risk loci to genes by their physical location using the Ensembl Variant Effect Predictor (VEP)³⁶. Of the lead variants, 59% (72/123) were within a transcript of a protein-coding gene, and 80% (99/123) of the loci contained at least one protein-coding gene within 20kb, and 93% (114/123) within 250kb (Supplementary Table 3). Of the 123 lead variants, 5 were missense variants (in genes PLCE1, MRGPRE, SERPINA1, ZBTB4 and ZNF462), and 40 more missense variants were in high LD (r²>0.6) with the lead variants (Supplementary Table 7a). Of note, three variants with a predicted high impact consequence on protein function were in high LD with the lead variants: (1) a stop gained variant (rs34358) with lead variant rs42854 (r²=0.85) in gene ANKDD1B, (2) a splice donor variant (rs66880209) with lead variant rs1472662 (r²=0.71) in RP11-420K8.1 and (3) a splice acceptor variant (rs11042902) with lead variant rs4910165 (r²=0.69) in MRVI1 (Supplementary Table 7b).

We used stratified LDSC (S-LDSC) to partition migraine heritability by 24 functional genomic annotations^37,38. We observed enrichment for ten categories (Supplementary Fig. 3 and Supplementary Table 8), with conserved regions showing the highest enrichment (11.2-fold; P=1.95×10⁻¹⁰), followed by coding regions (8.1-fold; P=1.36×10⁻³) and enhancers (4.2-fold; P=3.64×10⁻⁴).

Prioritization of candidate genes

We mapped the 123 lead variants to genes via expression quantitative trait locus (eQTL) association using GTEx v.8 (ref. ³⁹) and data repositories included in FUMA⁴⁰ at a false discovery rate (FDR) of 5% (Methods). The lead variants were cis-eQTLs for 589 genes (Supplementary Table 9), and variants in high LD with the lead variants were cis-eQTLs for an additional 624 genes (Supplementary Table 10). In total, 84% (103/123) of lead variants were cis-eQTLs for at least one gene. Tibial artery had the highest number (47/123) of lead variants as cis-eQTLs in GTEx v.8, and it was the only tissue type where the enrichment was statistically higher (P=6.37×10⁻⁶) than expected based on the overall number of cis-eQTLs per tissue reported by GTEx (Supplementary Fig. 4 and Supplementary Note).

To prioritize candidate genes for the risk loci, we applied two approaches based on GTEx v.8 expression data: fine-mapping of causal gene-sets by FOCUS⁴¹ (Supplementary Table 11a) and a transcriptome-wide association study (TWAS) by S-PrediXcan⁴² combined with colocalization analysis using COLOC⁴³ (Supplementary Table 11b).

With posterior probability (PP)>0.5, FOCUS found candidate genes for 82 loci and S-PrediXcan + COLOC supported colocalization for 52 loci (Supplementary Table 11c). In total, 73 genes in 46 loci were prioritized by both methods. MRC2 and PHACTR1 were the only genes that both methods prioritized with strong evidence (PP>0.99 for same tissue) and without any other gene prioritized within their loci.

Two of the new risk loci contain genes (CALCA/CALCB and HTR1F) whose protein products are closely related to targets of two migraine-specific drug therapies⁴⁴. We observe a convincing association at the chromosome 11 locus that contains the CALCA and CALCB genes encoding CGRP itself (lead SNP rs1003194, P=2.43×10⁻¹⁰; Fig. ​Fig.2a),2a), while none of the genes encoding CGRP receptor proteins (CALCRL, RAMP1 or RCP) show a statistically comparable association (all P>10⁻⁴; Supplementary Fig. 5). Variant rs1003194 is a cis-eQTL for CALCB, but also for COPB1, PDE3B and INSC (Supplementary Table 9) and FOCUS prioritizes CALCA, CALCB and INSC (Supplementary Table 11c). In addition, a new locus on chromosome 3 contains HTR1F (lead SNP rs6795209, P=1.23×10⁻⁸; Fig. ​Fig.2b),2b), which encodes the serotonin 5-HT_1F receptor. Variant rs6795209 is a significant cis-eQTL for HTR1F, as well as for three other genes (CGGBP1, ZNF654, C3orf38) in the same locus (Supplementary Table 9). FOCUS or S-PrediXcan + COLOC did not prioritize HTR1F based on gene expression data (Supplementary Table 11c).

Open in a separate window

Fig. 2

LocusZoom plots of two previously unknown migraine loci with genes that are targets of recent migraine-specific drugs.

a, Locus containing CALCA and CALCB genes, encoding CGRP, which is the target of preventive and acute therapies via monoclonal antibodies and gepants. b, Locus containing the HTR1F gene, which encodes a serotonin 5-HT_1F receptor that is the target of acute therapies via ditans. The xaxis shows the chromosomal location, and the yaxis shows the uncorrected two-sided negative log₁₀ (P) from the inverse–variance weighted fixed-effects meta-analysis with 102,084 cases and 771,257 controls. The squared correlation to the lead variant is shown by colors based on the UK Biobank data for variants that have an effective sample size ±20% of the lead variant’s effective sample size. Horizontal line corresponds to P=5×10⁻⁸. Blue graph shows the recombination rate.

Migraine subtypes with aura and without aura

Previously, Gormley et al.¹³ conducted subtype-specific GWAS with 6,332 MA cases against 144,883 controls and 8,348 MO cases against 139,622 controls, and reported that seven loci were GWS in MO but none were GWS in MA. Here, we added to the previous data 8,292 new MA and 6,707 new MO cases from headache specialist centers in Denmark and the Netherlands as well as from study collections in Iceland and UK Biobank (Table ​(Table2),2), for total sample sizes of 14,624 MA cases and 703,852 controls, and 15,055 MO cases and 682,301 controls. We estimated the effect size for each subtype at the 123 lead variants of the migraine GWAS (Supplementary Table 3b,c and Supplementary Data 4 and 5) and detected four GWS variants in the MA meta-analysis and 15 GWS variants in the MO meta-analysis. We also estimated a probability that the lead variant is either subtype-specific (that is, associated only with MO or with MA but not with both), shared by both subtypes, or not associated with either subtype (Methods; Supplementary Table 12a and Supplementary Data 6). With a probability above 95%, three lead variants (that is, rs12598836 in the HMOX2 locus, rs10405121 in the CACNA1A locus and rs11031122 in the MPPED2 locus) are MA-specific, while two lead variants (that is, rs7684253 in the locus near SPINK2 and rs8087942 in the locus near FECH) are MO-specific at a similar threshold. Nine lead variants were shared by MA and MO with >95% probability (Fig. ​(Fig.3a).3a). In addition to the five subtype-specific lead variants, four other lead variants also showed differences in effect size between the subtypes (P<0.05/123) (Fig. ​(Fig.3b3b).

Open in a separate window

Fig. 3

Lead variants stratified by migraine subtype for risk loci with MAF >5%.

a, Axes show the negative log₁₀ (P) of MO (xaxis; n=697,356, 15,055 MO cases and 682,301 controls) and MA (yaxis; n=718,476, 14,624 MA cases and 703,852 controls) analyses. Two-sided Pvalues are derived from inverse–variance weighted fixed-effect meta-analyses and are uncorrected. Symbols that are colored and annotated indicate >95% posterior probability that a nonzero effect is present in both MO and MA (model BOTH), or that the effect is present only in MO or only in MA but not both (models MO and MA, respectively). Variants with a probability <95% for each of the three models are shown as gray. Dashed lines show the genome-wide significance threshold (P=5×10⁻⁸). b, Axes show logarithm of odds ratios (OR) for MO (xaxis; n=697,356, 15,055 MO cases and 682,301 controls) and MA (yaxis; n=718,476, 14,624 MA cases and 703,852 controls) calculated for the migraine risk allele. The effects at variants that have been colored and annotated differ between the subtypes at significance level of 0.0004=0.05/123. The 95% CIs for the logarithm of odds ratios are shown for the annotated variants. Dashed lines show the coordinate axes and the diagonal.

Phenome-wide association scans with National Human Genome Research Institute GWAS Catalog and FinnGen R4

Next, we conducted phenome-wide association scans (PheWAS) for the lead variants for 4,314 traits with reported associations in the National Human Genome Research Institute (NHGRI) GWAS catalog (https://www.ebi.ac.uk/gwas/) and for the GWAS summary statistics of 2,263 disease traits in the FinnGen release 4 data. We identified 25 lead variants that were reported to be associated with 23 different phenotype categories (Methods) in the GWAS Catalog, and 17 lead variants with 26 defined disease categories in FinnGen at P<1×10⁻⁵. The categories with the highest number of reported associations were cardiovascular disease (7 lead variants) and blood pressure (6 lead variants) in the GWAS catalog, and diseases of the circulatory system (11 lead variants) in FinnGen. When we performed PheWAS for all variants in high LD (r²>0.6) with the lead variants, we observed associations for 79 loci with 54 different phenotype categories in the GWAS Catalog, and for 41 loci with 26 disease categories in FinnGen (Supplementary Table 13a and Supplementary Fig. 6).

These findings are consistent with previous results that migraine is a risk factor for several cardiovascular traits^45–47, and genetically correlated with blood pressure^48,49. However, we did not observe a trend in the direction of the allelic effects between migraine and coronary artery disease (CAD) or migraine and blood pressure traits (Supplementary Table 13d) using the latest meta-analysis of the CARDIoGRAMplusCD4 Consortium⁵⁰ (n=336,924) and blood pressure GWAS from UK Biobank⁵¹ (n=422,771).

Quality control

Before the meta-analysis, a standard quality control (QC) protocol was applied to each individual GWAS. Related individuals were removed from all other cohorts except HUNT (which modeled relatedness via a logistic mixed model) by using an identity by descent cut-off of 0.185 or smaller. Multiallelic variants were excluded from all studies, and only variants that satisfied the following thresholds were kept for further analysis: minor allele frequency (MAF)>0.01, IMPUTE2 info or MACH r²>0.6 and, when available, Hardy-Weinberg equilibrium (HWE) Pvalue >1×10⁻⁶ and missingness <0.05. Variants were matched by chromosome, position and alleles to the UK Biobank data. Indels were recoded as insertions (I) and deletions (D). For each study, SNPs with an effect allele frequency (EAF) discrepancy of >0.30 and indels with EAF discrepancy of >0.20 to UK Biobank were excluded. MAF and EAF plots of cohorts against the reference cohort are shown in Supplementary Data 7. We conducted a sensitivity analysis on strand-ambiguous SNPs (with alleles A/T or G/C), by counting, for each pair of studies, how often the same allele of A/T or G/C SNP was coded as the minor allele in both cohorts, as a function of MAF threshold (Supplementary Table 17). Minor alleles were same at least in 97.39% of the SNPs without MAF threshold and the corresponding proportions were 99.96% and 79.58% when MAF<0.25 and when MAF>0.4, respectively. The very high concordance for SNPs with MAF<0.25 suggests that the strand-ambiguous SNPs were labeled consistently for almost every SNP. Therefore, we did not exclude any SNPs based on possible labeling mismatches due to strand ambiguity.

Statistical analysis

All statistical tests conducted were two-sided unless otherwise indicated. The GWAS for the individual study cohorts were performed by logistic regression with an additive model of imputed dosage of the effect allele on the log-odds of migraine. The analyses for IHGC2016 (ref. ¹³) and 23andMe¹⁹ have been described before. For UKBB data and GeneRISK data, we used PLINK v.2.0 (ref. ⁸¹). For HUNT data, we used a logistic mixed model with the saddlepoint approximation as implemented in SAIGE v.0.20 (ref. ⁸²) that accounts for the genetic relatedness. All models were adjusted for sex and at least for the four leading principal components of the genetic population structure (Supplementary Table 18). Age was used as a covariate when available. A detailed description is provided in Supplementary Note. For the chromosome X meta-analysis, male genotypes were coded as {0,2} in all cohorts, and the GWAS were conducted with an X chromosome inactivation model that treats hemizygous males as equivalent to homozygous females⁸³.

We performed an inverse–variance weighted fixed-effect meta-analysis on the five study collections by using GWAMA⁸⁴. After the meta-analysis, we excluded the variants with effective sample size N_eff<5,000 to remove results with very low precision compared with most variants and were left with 10,843,197 variants surpassing the QC thresholds. We estimated the effective sample size for variant i as

where f_i is the effect allele frequency for variant i and s.e._i is the s.e. for variant i estimated by the GWAS software. This quantity approximates the value 2N t(1–t)I, where N is the total sample size (cases+controls), t is the proportion of cases and I is the imputation info (derivation in Supplementary Note).

Risk loci

There were 8,117 GWS variants with the meta-analysis Pvalue <5×10⁻⁸. For 8,067 of them that were available in UK Biobank, an LD matrix was obtained from UK Biobank using a random sample of 10,000 individuals included in the UKBB GWAS. We defined the index variants as the LD-independent GWS variants at LD threshold of r²<0.1 in the following way. First, the GWS variant with the lowest Pvalue was chosen and, subsequently, all GWS variants that were in LD with the chosen variant (r²>0.1) were excluded. Next, out of the remaining GWS variants, the variant with the lowest Pvalue was chosen and the GWS variants in LD with that variant were excluded. This procedure was repeated until there were no GWS variants left. Out of the 8,067 variants with LD information, 170 were LD-independent (at r²<0.1). For 18/50 variants that were not found in UK Biobank, LD information was available from the 23andMe data, and all 18 variants were in LD (r²>0.1) with some index variant. Of the 18 variants, 2 (rs111404218 and rs12149936) had lower Pvalue than the original index variant they were in LD with and, hence, they replaced the original index variants. For 32 GWS variants, LD remained unknown. Thus, at this stage, the GWS associations were represented by 202 = 168+2+32 index variants.

Next, to define the risk loci and their lead variants, an LD block around each index variant was formed by the interval spanning all GWS variants that were in high LD (r²>0.6) with the index variant. Sizes of these regions ranged from 1bp (only the variant itself, for example, the variants with unknown LD) to 1,089kb. Sets of regions that were less than 250kb away from each other were merged (distance from the end of the first region to the beginning of the second region). This definition resulted in 126 loci. All other GWS variants were included in their nearest locus based on their position and the locus boundaries were updated and, finally, loci within 250kb from each other were merged. This resulted in our final list of 123 risk loci. Each risk locus was represented by its lead variant defined as the variant with the lowest Pvalue and named by the nearest protein-coding gene to the lead variant or by the nearest noncoding gene if there was no protein-coding gene within 250kb. The term ‘Near’ was added to the locus name if the lead variant did not overlap with a gene transcript. We note that the nearest gene to the lead variant need not be a causal gene. None of the 32 variants without LD information became a lead variant of a risk locus because all had a variant in the vicinity with a smaller Pvalue.

We annotated and mapped these loci by their physical position to genes by using the Ensembl Variant Effect Predictor (VEP, GRCh37)³⁶. We used two different thresholds for annotating the nearest genes: a distance of 20kb and 250kb to the nearest transcript of a gene. The filtered results including all variants within a gene or a regulatory element are presented in Supplementary Table 7b.

Stepwise conditional analysis

We performed a stepwise conditional analysis (CA) on each risk locus by using FINEMAP v.1.4 (ref. ⁸⁵). FINEMAP uses GWAS summary statistics together with an LD reference panel and does not require individual-level data. When the reference LD does not accurately match the GWAS data, full fine-mapping is prone to false positives⁷⁹. A simpler stepwise CA is more robust to inaccuracy in reference LD because CA has a much smaller search space than full fine-mapping, and therefore CA is less likely to run into most problematic variant combinations where LD is very inaccurate. Since we did not have the full in-sample LD from our GWAS data, we carried out only the CA and not the full fine-mapping. For the CA, we included only the SNPs, but no indels, and we used the same reference LD from the UK Biobank data as we used to define the risk loci. We restricted the CA only to the variants with a similar effective sample size (N_eff) by using a threshold of ±10% of the N_eff of the lead SNP of the risk locus, because our summary statistics came from the meta-analysis where sample sizes per variant vary greatly. This filter excluded approximately 17% of all GWS variants and was necessary since otherwise CA led to spurious conditional Pvalues, such as P<10⁻²⁵⁰, for some loci. Consequently, for two of the loci where the lead variant was an indel, the lead variant was not included in the CA. For such regions, we checked that the new lead variant from the CA output was in LD (r²>0.3) with the original lead variant. For one locus (rs111404218) where the lead variant does not have LD information in the UK Biobank data, there were no GWS variants left in the CA after filtering by N_eff. We used the standard GWS (P<5 × 10⁻⁸) threshold to define the secondary variants that were conditionally independent from the lead variant. The CA results are in Supplementary Tables 6a,b.

eQTL mapping to genes and tissues

We used two data sources to map the risk variants to genes via eQTL associations. From the GTEx v.8 database (https://gtexportal.org), we downloaded the data of 49 tissues. We first mapped all 123 lead variants to all significant cis-eQTLs across tissues using the FDR cut-off of 5% as provided by the GTEx project³⁹. Next, we also mapped the variants in high LD (r²>0.6) with the lead variants to all significant cis-eQTLs. Finally, we filtered the results to include only the new significant gene-tissue pairs that were not implicated by the lead variants. Results are shown in Supplementary Tables 9 and 10.

With FUMA v.1.3.6 (ref. ⁴⁰), we mapped the 123 lead variants, and the variants in high LD (r²>0.6) with the lead variants, to the other eQTL data repositories provided by FUMA except GTEx, that is, Blood eQTL Browser⁸⁶, BIOS QTL browser⁸⁷, BRAINEAC⁸⁸, MuTHER⁸⁹, xQTLServer⁹⁰, CommonMind Consortium⁹¹, eQTLGen⁹², eQTL Cataloque⁹³, DICE⁹⁴, scRNA eQTLs⁹⁵ and PsychENCODE⁹⁶. Results are shown in Supplementary Tables 9 and 10.

To study whether the lead variants were enriched in any of the 49 tissues from GTEx v.8, we fitted a linear regression model where the number of lead variants that are significant cis-eQTLs for a specific tissue was used as the outcome, and the overall number of genes with at least one significant cis-eQTL reported by GTEx for the tissue was the predictor³⁹. We did a separate regression model for each tissue type by leaving the tissue of interest out from the model, and we used the model fitted on the other tissues for predicting the outcome variable for the tissue type of interest. Finally, we checked in which tissues the true observed number of migraine lead variants was outside of the 95% prediction intervals as given by the function ‘predict.lm(‘interval=’prediction’)’ in R software. Details of the procedure are in the Supplementary Note.

LD-score regression

We estimated both the SNP heritability ($h_{\mathrm{SNP}}^2$) of migraine and pairwise genetic correlations (r_G) between each pair of study collections using LDSC v.1.0.0 (refs. ^30,35). SNP heritability and genetic correlations were estimated using European LD scores from the 1000 Genomes Project Phase3 data for the HapMap3 SNPs, downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/. We reformatted the meta-analysis association statistics to LDSC format with munge-tool, which excluded variants that did not match with the HapMap3 SNPs, had strand ambiguity (that is, A/T or G/C SNPs), MAF<0.01 or missingness more than two-thirds of the 90th percentile of the total sample size, or resided in long-range LD regions⁹⁷, in centromere regions or in the major histocompatibility locus (MHC) of chromosome 6, leaving 1,165,201 SNPs for the LDSC analyses. We used a migraine population prevalence of 16% and a sample proportion of cases of 11.7% = 102,084/(102,084+771,257) to turn the LDSC slope into the estimate of $h_{\mathrm{SNP}}^2$ on the liability scale⁹⁸. Pairwise genetic correlation results are listed in Supplementary Table 2. We note that in the previous migraine meta-analysis¹³, LDSC reported $h_{\mathrm{SNP}}^2$ value of 14.6% (13.8–15.5%), which was considerably larger than the value 11.2% (10.8–11.6%) that we report in our analysis. When we ran our LDSC pipeline on the data of Gormley et al.¹³, we estimated $h_{\mathrm{SNP}}^2$ value of 10.6% (10.1–11.1%). Thus, it seems that our liability transformation estimates lower values of heritability than the transformation used by Gormley et al.¹³.

Stratified LD-score regression

We used S-LDSC to partition the SNP heritability by functional genomic annotations³⁷. We used the baseline-LD model³⁸ that contains 75 annotations, including conserved, coding and regulatory regions of the genome and different histone modifications. The baseline-LD model adjusts for MAF- and LD-related annotations, such as recombination rate and predicted allele age, which decreases the risk of model misspecification^37,38,99. We used the same QC as with the univariate LDSC, and the baseline LDv.1.1 European LD scores estimated from the 1000 Genomes Project Phase3, downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/. We set the significance threshold for enrichment of individual binary functional annotations to α=0.05/24, as we considered only 24 unique functional annotations without the flanking regions. Results are listed in Supplementary Table 8.

Subtype analyses of MA and MO

First, we combined new MA and MO data (Table ​(Table2)2) with the previously used migraine subtype-specific meta-analysis data¹³, and estimated migraine subtype-specific effect sizes for the 123 lead variants from the migraine meta-analysis. We tested how often the direction of allelic effects was similar between the IHGC MA/MO and the new cohorts using a binomial test (Supplementary Table 12b). Next, we stratified the lead variants by using the information from the migraine subtype-specific analyses. For each of the variants, we estimated probabilities between four possible explanations of the observed data that we call ‘NULL’, ‘MO’, ‘MA’ and ‘BOTH’. Under model NULL, the effect is not present in either of the migraine subtypes (that is, the effect is zero); under model MO or MA, the effect is present only in MO or only in MA but not in both; and under model BOTH, a nonzero effect is shared by both MO and MA. We used a Bayesian approach for model comparison that combines a bivariate Gaussian prior distribution on the two effect sizes with a bivariate Gaussian approximation to the likelihood using GWAS summary statistics¹⁰⁰. Across all models, the prior s.d. for the effect is 0.2 on the log-odds scale for nonzero effects and 0 for a zero effect. The bivariate priors for the four models are as follows: NULL assumes a zero effect in both migraine subtypes, MO and MA assume a nonzero effect for one subtype and a zero effect for the other subtype, and BOTH combines the fixed-effect model (exactly the same effect in both subtypes) with the independent-effects model (the two effect sizes are nonzero but uncorrelated with each other) with equal weights. Finally, we assumed that each of the four models (NULL, MO, MA and BOTH) is equally probable a priori, which we considered an appropriate assumption since all these variants show a convincing association to overall migraine (P<5×10⁻⁸). Then we used the Bayes formula to work out the posterior probability on each model. The results are shown in Fig. ​Fig.3a,3a, thresholded by a probability cut-off of 95%, and in Supplementary Table 12a. The correlation parameter between MO and MA GWAS statistics needed in the bivariate likelihood approximation was estimated to be 0.148 using the empirical Pearson correlation of the effect size estimates of the common variants (MAF>0.05), which did not show a strong association to either of the migraine subtypes (P>1×10⁻⁴)¹⁰¹. We tested whether the effect sizes between MA and MO were equal at a Bonferroni corrected significance threshold of α = 0.05/123 by using a normal approximation and accounting for the correlation in effect size estimators.

We note that the amount of information in the data (‘statistical power’) is taken into account automatically in this model comparison, which we consider an advantage compared with a comparison of the raw Pvalues between the subtype analyses that does not automatically account for statistical power. In particular, observing a GWS Pvalue (P<5×10⁻⁸) in one subtype but not in the other subtype is not yet evidence for a subtype-specific locus, because the effect could still be nonzero also for the other subtype but simply lack power to reach the stringent GWS threshold. Finally, we point out that the inference in the model comparison approach is conditional on the particular set of models being included in the comparison, as well as on the particular choice of the prior distributions.

PheWAS with NHGRI GWAS catalog and FinnGen R4

We performed PheWAS for the 123 lead variants using the NHGRI GWAS catalog and the FinnGen R4 GWAS summary statistics. In addition, we performed the same lookups for the 123 risk loci including all variants in high LD (r²>0.6) with the lead variants. With the GWAS catalog, we first downloaded all the available results (4,314 traits) from the GWAS catalog webpage (accessed April 6, 2020). Next, we obtained all the associations for the 123 risk loci with all the high LD variants included using Pvalue thresholds of P<1 × 10⁻⁵, P<1 × 10⁻⁶ and P<1 × 10⁻⁴ (Supplementary Table 13a–c). Because the GWAS catalog includes results from several different GWAS for the same phenotype or for a very similar phenotype with a different name, we divided the phenotype associations into broader categories. The new categories are listed in Supplementary Table 19. The same approach was used for the PheWAS of FinnGen R4. We first downloaded all the available summary statistics (2,263 endpoints) and, next, obtained all the associations for the 123 risk loci using the same three Pvalue thresholds as with the GWAS catalog (Supplementary Table 13a–c). We also divided similar endpoints into broader categories, which are listed in Supplementary Table 20.

We tested the direction of allelic effects between migraine and the following three traits that shared multiple associated variants with migraine: CAD⁵⁰, diastolic blood pressure⁵¹ and systolic blood pressure⁵¹. We first took all migraine lead variants that were available also in the summary statistics of the other trait without any Pvalue threshold and used a binomial test to test whether the proportion of variants with same direction of effects was 0.5. Next, we used a Pvalue threshold of 1×10⁻⁵ for the association with the other trait. Results are in Supplementary Table 13d.

LD-score regression applied to specifically expressed genes

We used LDSC-SEG¹⁴ to identify tissues and cell types implicated by the migraine GWAS results. LDSC-SEG uses gene expression data and GWAS results from all variants together with an LD reference panel. For our analyses, we used the same QC as for the other LDSC analyses and six different sets of readily constructed annotation-specific LD scores downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/LDSC_SEG_ldscores/: multitissue gene expression, multitissue chromatin, GTEx brain, Cahoy, Corces ATAC and ImmGen LD Scores. FDR was controlled by the Benjamini–Hochberg method. The results are in Supplementary Table 14a–f. There were no significant results with the Cahoy, Corces ATAC and ImmGen data at FDR 5%.

Multimarker analysis of genomic annotation

We applied multimarker analysis of genomic annotation (MAGMA) v.1.09 (ref. ⁵³) to identify genes and gene sets associated with the migraine meta-analysis results. First, we mapped the meta-analysis SNPs to 18,985 protein-coding genes based on their physical position in the National Center for Biotechnology Information 37 build by using default settings of MAGMA. Next, we performed a gene-based analysis using the default SNPwise-mean model and the same UK Biobank LD reference as for the other analyses. We applied a Bonferroni correction (α=0.05/18,985) to identify significantly associated genes for migraine with the results listed in Supplementary Table 16a. Finally, we used the results from the gene-based analysis to perform a gene-set analysis by using two different gene-set collections from the Molecular Signature Database v.7.0 (refs. ^102,103): the curated gene sets containing 5,500 gene sets and the GO gene sets containing 9,988 gene sets. We performed the gene-set analysis using the competitive gene-set model and one-sided test that tests whether the genes in the gene-set are associated more strongly with the phenotype compared to the other genes. To correct for multiple testing, we used a Bonferroni correction (α=0.05/(5,500+9,988)). Results are in Supplementary Table 16b,c and in Supplementary Fig. 7.

DEPICT

DEPICT⁵² is an integrative tool to identify the most likely causal genes at associated loci, and enriched pathways and tissues or cell types in which the genes from the associated loci are highly expressed. As an input, DEPICT takes a set of trait-associated SNPs. First, DEPICT uses coregulation data from 77,840 microarrays to predict biological functions of genes and to construct 14,461 reconstituted gene sets. Next, information of similar predicted gene functions is used to identify and prioritize gene sets that are enriched for genes in the associated loci. For the tissue- and cell-type-enrichment analysis, DEPICT uses a set of 37,427 human gene expression microarrays. We used DEPICT v.1.194 and ran the analyses twice for each of the Pvalue thresholds for clumping, as recommended⁵², and using the default settings of 500 permutations for bias adjustment and 50 replications for the FDR estimation and for the Pvalue calculation. As an input, we used only the autosomal SNPs and the same UK Biobank LD reference data as for the other analyses. First, we ran the analysis using a clumping Pvalue threshold of 5×10⁻⁸ that resulted in 165 clumps formed from 7,672 variants (Supplementary Table 15d–f). Second, we used a Pvalue threshold of 1×10⁻⁵ leading to 612 clumps formed from 22,480 variants (Supplementary Table 15a–c).

Transcriptome-wide association study and colocalization

We performed a transcriptome-wide association study (TWAS) by S-PrediXcan⁴² v.0.7.5 using GTEx v.8 multivariate adaptive shrinkage models (MASHR-M) for 49 tissues downloaded from predictdb.org and the European 1000 Genomes v.3 LD reference panel (hg38; https://zenodo.org/record/3657902/). We followed the recommended QC protocol, and first harmonized and imputed the migraine summary statistics to ensure an optimal overlap with the GTEx v.8 expression weights. After harmonization and summary statistic imputation, 8,909,736 variants were available for the TWAS. We performed the analysis with default settings to identify significant gene-tissue pairs. We applied a Bonferroni corrected significance level of α=0.05/662,726, corresponding to the number of unique gene-tissue pairs tested.

Next, we performed colocalization analysis with COLOCv.4.0.4(ref. ⁴³) R package for the 1,844 significant gene-tissue pairs to indicate pairs that could be due to LD contamination. COLOC compares five hypotheses where the null hypothesis (H0) corresponds to no association to either eQTL or GWAS, H1 and H2 correspond to associations with only one of the traits, H3 corresponds to association with both eQTL and GWAS but at distinct causal variants, and H4 corresponds to association with both eQTL and GWAS at a shared causal variant. We set a prior probability for colocalization as p₁₂=5×10⁻⁶ for all tested regions and restricted the analysis to variants that had N_eff±10% of the N_eff of the lead variant of the region. Results are presented in Supplementary Table 11b.

Fine-mapping of causal gene sets

To prioritize genes for the migraine loci, we applied a gene-based fine-mapping approach using fine-mapping of causal gene sets (FOCUS) v.0.7 (ref. ⁴¹). FOCUS is a Bayesian approach that models predicted expression correlations among TWAS signals to estimate posterior probabilities for all genes within a tested region.

We used the European 1000 Genomes v.3 LD reference panel and same GTEx v.8 predicted expression weights for the 49 tissues as with S-PrediXcan. First, we mapped the migraine summary statistics from hg37 to hg38 with UCSC liftOver¹⁰⁴. Next, we followed the suggested QC protocol and applied the modified munge-tool to obtain cleaned summary statistics. After the QC steps, we had 6,237,177 variants left for the analysis. We performed tissue-prioritized fine-mapping of gene-sets for the 49 tissues with otherwise default settings except that we increased the Pvalue threshold to 1×10⁻⁴ so that the fine-mapping would cover most of the same regions that contained at least one significant gene-tissue pair by S-PrediXcan. Posterior inclusion probability (PIP) from FOCUS is reported for all available significant S-PrediXcan gene-tissue pairs in Supplementary Table 11b, and all prioritized genes by FOCUS with PIP>0.9 are reported in Supplementary Table 11a.

We thank the study participants for their contribution to this research. We also thank the numerous individuals who contributed to sample collection, storage, handling, phenotyping and genotyping for each of the individual cohorts. We acknowledge the participants and investigators of the FinnGen study. This research has been conducted using the UK Biobank Resource under Application Number 22627. We are supported by following grants: the US National Institute of Neurological Disorders and Stroke (NINDS) of the US National Institutes of Health (NIH) (grant numbers R21NS09296 and R21NS104398 (D.I.C.)), the Finnish innovation fund Sitra and Finska Läkaresällskapet (E.W.), the Academy of Finland (grant nos. 288509, 312076, 336825 (M.P.)), the Sigrid Juselius Foundation (M.P. and S.R.), the Academy of Finland Center of Excellence in Complex Disease Genetics (grant no. 312062 (S.R.)), the Finnish Foundation for Cardiovascular Research (S.R.), University of Helsinki HiLIFE Fellow and Grand Challenge grants (S.R.), The Novo Nordisk Foundation (NNF14CC0001 and NNF17OC0027594 (T.F.H. and K.B.)), CANDY foundation (CEHEAD) (T.F.H.), and the South-Eastern Norway Regional Health Authority (grant no. 2020034 (B.S.W.)). A list of study-specific acknowledgements and funding information can be found in the Supplementary Note.