GWAS in childhood acute lymphoblastic leukemia reveals novel genetic associations at chromosomes 17q12 and 8q24.21

Childhood acute lymphoblastic leukemia (ALL) (age 0–14 years) is 20% more common in Latino Americans than non-Latino whites. We conduct a genome-wide association study in a large sample of 3263 Californian children with ALL (including 1949 of Latino heritage) and 3506 controls matched on month and year of birth, sex, and ethnicity, and an additional 12,471 controls from the Kaiser Resource for Genetic Epidemiology Research on Aging Cohort. Replication of the strongest genetic associations is performed in two independent datasets from the Children’s Oncology Group and the California Childhood Leukemia Study. Here we identify new risk loci on 17q12 near IKZF3/ZPBP2/GSDMB/ORMDL3, a locus encompassing a transcription factor important for lymphocyte development (IKZF3), and at an 8q24 region known for structural contacts with the MYC oncogene. These new risk loci may impact gene expression via local (four 17q12 genes) or long-range (8q24) interactions, affecting function of well-characterized hematopoietic and growth-regulation pathways.

A cute lymphoblastic leukemia (ALL), the most common type of cancer in children, has an increasingly clear genetic etiology with a role for both rare high-penetrance genetic predisposition alleles 1 and more common low-penetrance alleles. These latter alleles were discovered via genome-wide association studies (GWAS), which identified genetic risk factors at ARID5B, CEBPE, IKZF1, CDKN2A, PIP4K2A, and GATA3 2-5 . Additional loci were recently identified at LHPP and ELK3 6 . While most of these genes encode hematopoietic transcription factors (TFs), the causal variants underlying these genetic associations remain unclear, apart from a missense polymorphism in CDKN2A 7-9 and single-nucleotide polymorphisms (SNPs) disrupting TF-binding sites in the promoters of CEBPE 10 and ARID5B 11 . These prior GWAS focused on non-Latino white populations; however, the highest rate of childhood ALL globally is found in Latinos in the United States (US) 12 . We sought to identify additional risk alleles involved in the etiology of ALL in both Latinos and non-Latinos by leveraging a large, ethnically diverse population in California in conjunction with an array specifically designed for the Latino population, while supplementing with additional accessible controls genotyped on the same platform. We replicated our findings in two independent childhood ALL case-control datasets from the Children's Oncology Group (COG) and the California Childhood Leukemia Study (CCLS). We identify two new genetic loci impacting childhood ALL risk, which will aid in our understanding of the etiology of this disease.

Results
Association analysis. After quality control, 757,935 polymorphic autosomal SNPs and 19,240 subjects (3263 cases, 15,977 controls, which included genotyped controls run specifically for this study and additional public controls, see Supplementary Table 1) remained for discovery analyses. Case-control association tests were stratified by race/ethnicity including 10,533 Latino subjects (1949 cases, 8584 controls), 4735 non-Latino white subjects (1184 cases, 3551 controls), and 3972 African-American subjects (130 cases, 3842 controls), and then the race/ethnicity-stratified GWAS were combined via meta-analysis. Sample size and ancestral heterogeneity precluded the inclusion of more minor populations (e.g., Southeast Asians, Indian subcontinental, and other races/ethnicities) for the current analysis. Associations reached genome-wide statistical significance (i.e., P < 5.0 × 10 −8 in an additive logistic regression analysis) at known ALL risk loci, including: ARID5B, IKZF1, PIP4K2A, and CEBPE among Latinos, and ARID5B, IKZF1, and PIP4K2A for non-Latino whites (Supplementary Table 2). No SNPs reached genome-wide significance among African Americans, likely due to the relatively small sample size. Analyses by each racial/ethnic subgroup did not reveal any novel associations at genome-wide statistical significance.

PIP4K2A
CDKN2A CEBPE IKZF3 SP4 8q24 Fig. 1 Manhattan plot of P values resulting from meta-analysis of three separate analyses in the California population (Latinos, non-Latino whites, and non-Latino blacks) including study-specific controls and GERA controls (total N = 3263 cases, 15,977 controls). The locations of prior identified GWAS hits are indicated by their gene acronyms. New genetic associations identified in this study are highlighted in bold. These association peaks exceed the Bonferroni genome-wide P value cutoff of 5 × 10 −8 indicated by a horizontal line Latino controls from the CCLS (a set of cases and controls from California distinct from our discovery set). Using these data, we successfully replicated the novel SNP associations at 8q24 (OR = 1.30, 95% CI: 1.14-1.49, P meta = 1.29 × 10 −4 ) and near IKZF3 (OR = 1.14, 95% CI: 1.03-1.26, P meta = 0.013) ( Table 1). For some datasets, the top hit SNP was imputed as shown (Table 1). While the SP4 locus did not replicate in the combined COG European and CCLS Latinos replication datasets (P meta = 0.075), it was nominally associated with ALL risk in the COG sample, with a similar direction and magnitude of effect (OR = 1.15, 95% CI: 1.02-1.31, P = 0.025). SNP imputation enabled improved resolution of the associated regions in discovery samples (Fig. 2a-c). The SP4 association peak is confined to the SP4 gene only (Fig. 2a) and conditional analyses adjusting for rs2390536 revealed no independent risk alleles (Supplementary Figure 1). Because rs2390536 did not replicate in our combined replication set of CCLS Latinos and COG Europeans, it may be a false-positive signal in our discovery GWAS and therefore is not explored further here.
The chromosome 8q24 association peak was located within a 2.5 Mb region dense with genome-wide significant loci for a variety of developmental traits, including: cleft lip 14 , hypospadias 15 , and infant length 16 , as well as several cancers 17 . However, our ALL association peak was located within a 100 kb subregion (Fig. 2b), which has not been previously associated with any other phenotype. Conditional analyses adjusting for the lead SNP in the region greatly attenuated SNP associations in the region (Supplementary Figure 2), with the most significant residual association signal found at rs35380634, located some 300 kb telomerically (P conditional = 3.29 × 10 −5 ).
The ALL association peak on chromosome 17q12 is broad, covering approximately 200 kb and encompassing at least six genes, including IKZF3 (Fig. 2c). The lead SNP, rs2290400, is also a known GWAS hit for type 1 diabetes and asthma 18,19 , although the ALL risk allele is protective against development of these other immune-related phenotypes. Conditional analyses adjusting for rs2290400 effectively eliminated any residual association in the region, indicating that additional risk loci proximal to the primary peak likely do not exist (Supplementary Figure 3). Because IKZF1 also contains ALL risk alleles identified in previous GWAS, we tested for the presence of statistical interaction between our top genotyped SNP in IKZF1 (rs11978267; P CCRLP = 2.15 × 10 −29 ) and our top genotyped SNP in IKZF3. We detected no significant epistatic effect in either case-control analyses (P = 0.62) or case-only analyses (P = 0.16).
Our study is designed to discover SNPs associated with case status rather than the clinical course of disease. Within the COG replication dataset, information on patient relapse was available for all 959 cases, allowing for an assessment of risk of relapse. One hundred forty-three patients experienced a relapse, while 816 did not. Logistic regression analyses among these cases did not suggest that either the rs4617118 risk allele on 8q24.1, nor the rs2290400 risk allele near IKZF3, was associated with an increased risk of relapse (OR = 1.03 and OR = 0.89; P = 0.87 and P = 0.39, respectively).
Functional assessments. Top hit SNPs in the 8q24 and IKZF3 region, along with those that displayed larger effect sizes and P values within one order of magnitude, were annotated for functional evidence using HaploReg, RegulomeDB, and GTEx (Supplementary Tables 4 and 5). ENCODE-described motifs that bind at SNP locations, RegulomeDB scores, and other GWAS associations are indicated in the tables. Additional expression quantitative trait loci (eQTL) analysis was carried out using RNA sequencing (RNA-seq) data from lymphoblastoid cell lines (LCLs) in the gEUVADIS dataset (www.geuvadis.org).
8q24: The risk locus at 8q24 is located within a gene desert and is most proximal to, but outside of, a noncoding RNA LINC00977 (Fig. 3). The lead SNP, rs4617118, and those SNPs in tight linkage were not eQTL for any proximal or distal gene in any tissue in GTEx, nor for lymphoblastoid cells in gEUVADIS data. Many GWAS association peaks for other cancers and immune conditions exist on chromosome 8q24. Using chromatin conformation assays on LCLs 20 , the region exhibits two large structural chromosomal domains with MYC in the center (Fig. 3). SNP rs4617118 exists within one of these domains. This domain shows evidence of chromosome looping between MYC and the region containing rs4617118; this region also contains the gene GSDMC, on the distal side of rs4617118 (Fig. 3). GSDMC is a gasdermin family member, another of which (GSDMB) is also present within the chromosome 17q12 peak as explained below.
17q12: The genetic association peak at 17q12 was broad. The top SNP rs2290400 is a highly significant eQTL (P values < 10 −9 ) for two genes, GSDMB and ORMDL3, in whole blood (as found on gTEX, Fig. 4). An additional eQTL was identified at ZPBP2 in the LCL RNA-seq data (−log 10 (P value) of 5.737). Structural Hi-C analysis (using LCL data) did not reveal obvious structural interactions. The variant did not display any evidence of altering expression of IKZF3 using gTEX or gEUVADIS resources, though the linkage peak clearly covers this gene. Instead, the rs2290400 minor allele, protective against ALL, was clearly associated with decreased expression of its other associated gene targets GSDMB and ORMDL3. Further analysis revealed that rs2290400 disrupts a PRDM1 (BLIMP1) binding motif (Supplementary Table S4).

Discussion
Using a population-based case-series with both matched controls and additional controls from the same source population (albeit older in age) and two replication datasets, we identified and replicated two new genome-wide significant risk alleles for childhood ALL. These alleles augment our current understanding of the genetic etiology of childhood ALL by (i) adding one haplotype that is a strong eQTL for three genes and located at a wellestablished hematopoietic developmental gene (IKZF3), and (ii) adding an additional locus in a gene-poor region of chromosome 8q24 that is rich with genetic variants affecting risk for myriad different traits. A third association that reached genome-wide significance in our discovery set (at the gene SP4) did not replicate in the combined COG and CCLS replication data (P = 0.075).
Although it is not a confirmed risk factor for childhood ALL, it did reach an overall P value below the commonly accepted genome-wide significance level of 5 × 10 −8 in our meta-analysis and should be assessed in additional datasets. While the associated haplotype on 17q12 clearly contains IKZF3, the lead SNP rs2290400 was not identified as an eQTL for IKZF3 itself using available databases. Instead, these databases identified the neighboring genes ORMDL3, GSDMB, and ZPBP2 as top eQTLs. It is possible that the haplotype imparts more than a single genic effect on risk. Prior GWAS have shown that the ALL risk allele in rs2290400 (T) is protective against asthma 19 and type 1 diabetes 18 , suggesting a role for this SNP in immune regulation. An immunologic response role for the 17q12 haplotype is further supported by its association with a variety of other conditions of aberrant immune function, including: primary biliary cirrhosis 21 , systemic lupus erythematosus 22 , Crohn's disease 23 , and ulcerative colitis 24 .
It is intriguing that the well-characterized ALL risk allele in IKZF1 (Ikaros) is also protective against type 1 diabetes 25 , paralleling our findings at rs2290400 near IKZF3. The seeming paradox wherein a risk allele for asthma and diabetes is associated with decreased risk of childhood ALL, a disease also linked to immune activation (as exemplified by an association with documented history of infectious diseases) 26,27 , may be explained by complex interactions. While medical diagnosis of infection is associated with an increased risk of childhood ALL, exposure to infectious agents via daycare contact is associated with decreased risk 28 . However, it is important to note that IKZF3, ZPB2, ORMDL3, GSDMB, and GSDMA are structured in a large haplotype with co-regulatory features 29,30 . The fact that IKZF3 encodes the TF Aiolos, closely related to another ALL-associated TF Ikaros (IKZF1) that is also critical in hematopoietic development and cell fate, and that IKZF3 has been identified as a recurrent somatic deletion in leukemia 31 , suggests that it may be the functional target of our GWAS association here. In support of a causal role for rs2290400, this variant disrupts a binding motif for the B lymphocyte maturation-induced protein 1 (BLIMP1), a TF and known repressor of T-cell activation 32 . It is important to note that BLIMP1 is expressed in early stages of B-cell development 32 , and the SNP may affect risk via aberrant expression control of IKZF3 at critical hematopoietic differentiation junctures, similar to the activities of other SNPs known to affect B-cell leukemia risk 10 . If the ALL risk allele at rs2290400 favored B-cell differentiation while the alternate allele favored Tcell differentiation, this could potentially account for the inverse association with ALL vs. asthma and type 1 diabetes. We cannot at this point ascertain whether the 17q12 locus exerts functional effects via its impact on hematopoietic differentiation or immune response or both, and hence further research will be required. The risk locus at 8q24, centered on rs4617118, is a unique GWAS hit in a region rich with genetic associations. It is located 210,000 bases downstream of rs987525, a GWAS SNP for cleft palate 14 , and over 400,000 bases upstream of a cluster of SNPs associated with hypospadias, glioma, and white blood cell counts (GWAS catalog). The 8q24 SNP associated with white blood cell counts (rs1991866) is located near GSDMC. A related gene on chromosome 17 (GSDMB) is found within the leukemiaassociated haplotype block containing IKZF3. These genes are part of a family of genes involved in pyroptosis, or inflammatory programmed necrosis, which is a recently identified pathway involved in cell proliferation and differentiation, inflammatory remodeling, and cancer-related pathways 33 . Whether any relationship exists between the 8q24 locus and these pathways is unknown. Our preliminary suggestions on function of this allele center around a putative impact on MYC regulation, in concert with many other SNPs within the gene desert surrounding both sides of MYC which appear to disrupt long-range enhancer activity 34 , but do not exclude a potential effect on the function of GSDMC, which is harbored within the same looping domain.
Limitations of the current study include a lack of precise immunophenotype and tumor genetic information in our discovery population. Also, the use of imputed SNPs for some replication analyses as noted in Table 1 may introduce error, although we note that the information scores of our top hit SNPs were high. In addition, our initial goal was to discover aspects of the genetic etiology that explained the higher rate of childhood ALL in Latinos compared to non-Latino whites and the two new loci discovered here do not appreciably account for differences in genetic risk between these populations. Strengths of our discovery population include a population-based sampling strategy for both cases and controls that was not prone to participation bias, and the close match between cases and controls in terms of birth date, sex, and race/ethnicity are evidence for the lack of bias. Our results suggest that genetic risk factors of reasonable strength may continue to be discovered in population-based samples of childhood ALL.

Methods
Study subjects. The study protocol was approved by the Institutional Review Boards at the California Health and Human Services Agency, University of California (San Francisco and Berkeley), and Yale University. The California Department of Public Health, Genetic Diseases Screening Branch obtains blood samples (heel-prick blood spots) from all neonates born within the state for the purpose of genetic screening. The blood spots that remain after genetic screening have been archived at −20°C since 1982 and made available for appropriate research. We linked statewide birth records maintained by the California Department of Public Health (for the years of 1982-2009) to cancer diagnosis data from the California Cancer Registry (CCR, for the years of 1988-2011). Beginning in January 1988, California law has required all new cancer cases diagnosed in state residents to be reported to the CCR. Included in this analysis were children born in California during 1982-2009 and diagnosed with ALL at the age of 0-14 years per CCR record. Children who were born in California during the same period and not reported to CCR as having any childhood cancer were considered potential controls. For each case of childhood ALL, up to four control subjects were randomly selected from the pool of potential controls and matched to the case on year and month of birth, sex, and race/ethnicity (non-Latino white, non-Latino black, Latino, Asian/Pacific Islander, other). This database of cancer cases and controls is defined as the California Cancer Records Linkage Project (CCRLP). For the current GWAS, a newborn blood spot was pulled from the archive for each case, as well as one control randomly chosen from up to four matched controls for that case. Therefore, active genotyping for CCRLP specifically for this project involved an approximately equal number of cases and controls (see Supplementary Table 1).
CCRLP sample preparation and genotyping. DNA was extracted from neonatal blood spots according to the Qiagen DNA Investigator blood card extraction protocol. Three 3 mm punches (Wallac Multipuncher) were cut for each subject and subjected to automated DNA extraction (Qiacube). DNA specimens were assigned randomly to genotyping plates (93 unique samples per plate) and checked for randomization to ensure that case-control status, ethnicity, and sex were distributed evenly. DNA was genotyped on the Affymetrix Axiom World Array (Latino, or LAT) (plates 1-51, 60% of subjects) array 35 and the same array supplemented with MHC and KIR content for immunogenomic imputation purposes 36 on remaining samples (plates 52-86, 40% of subjects). The current manuscript contains analyses of array content common to both array designs, totaling 813,036 variants (i.e., the LAT array base content), plus an additional 8951 variants typed exclusively on plates 52-86. DNA samples were genotyped on an Affymetrix TITAN system, and raw image files were processed with Affymetrix Genetools to call genotypes. Plates 1-51 were clustered separately from plates 52 to 86 due to the differing SNP content, and then called genotypes were merged. Genome freeze hg19 (NCBI 37) was used throughout as the reference.
CCRLP genotyping quality control. In addition to randomizing samples, we included two duplicate samples per plate (174 total), with average genotype concordance >99%. SNP filtering was performed separately for the two array designs (i.e., plates 1-51 underwent SNP filtering separately from plates 52 to 86. Call-rate filtering for SNPs and samples were performed iteratively, as follows: SNPs with call rates <92% were removed, then samples with call rates <95% were excluded, then SNPs with call rates <97% were removed, then samples with call rates <96% were excluded. Additionally, 15 CEPH trios (parents and child) were genotyped to detect and exclude SNPs showing >2 Mendelian transmission errors. Any SNP that failed a quality control check on plates 1-51 was also excluded from analysis of plates 52-86, and vice versa. Within the merged dataset of all CCRLP cases and controls, we excluded any SNP displaying significant departure from Hardy-Weinberg equilibrium (P < 1.0 × 10 −4 ) among European-ancestry controls.
In addition to excluding samples due to call-rate filters (N = 167), samples with mismatched reported vs. genotyped sex were also excluded (N = 26). We performed identity-by-descent analyses in Plink on the merged dataset of all CCRLP cases and   37 . We incorporated genome-wide SNP array data from 1184 HapMap Phase 3 samples and performed multidimensional scaling analysis using unlinked autosomal biallelic SNPs with MAF > 0.05 allowing us to group 89.4% of cases and controls into one of three analytic groups: Latinos, non-Latino whites, and African Americans. Using hybridization intensity (log R ratio) data from array SNPs on chromosome 21, we identified 109 subjects with Trisomy 21 (~2.9% of total cases) and excluded these subjects from association tests due to the known differences in leukemogenesis in diploid vs. constitutively aneuploid patients 38 . Three additional cases were excluded because, on further examination, it was determined that ALL was not their primary cancer diagnosis.
Incorporation of additional public controls. Axiom genotyping data for a total of 73,298 non-Latino white, Latino, and African-American subjects from the Kaiser Resource for Genetic Epidemiology Research on Aging (GERA) Cohort were downloaded from dbGaP (Study Accession: phs000788.v1.p2). All GERA participants are members of Kaiser Northern California Health Care, similar to the CCRLP population. We restricted analyses to only those samples prepared using the Axiom "Type O" Reagent Kit (N = 12,641), as this was the same kit used for CCRLP samples. We excluded any SNPs that were not retained in the cleaned CCRLP dataset, any SNP with an MAF <5%, any SNP associated with control source (CCRLP vs. GERA) at P < 5 × 10 −8 in race/ethnicity-stratified analyses, and any SNP displaying significant departure from Hardy-Weinberg equilibrium (P < 1.0 × 10 −4 ) among European-ancestry controls. Samples with call rates <97% were also removed from analyses. The GERA control subject data were then merged with the CCRLP dataset and IBD analyses were performed in Plink, excluding one member of any sample pair that had an IBD proportion >0.20 (N = 161).
Discovery association analyses. Following SNP filtering procedures, a total of 757,935 polymorphic autosomal SNPs remained for analyses. Following sample filtering procedures, a total of 3263 cases (1822 males, 1441 females) and 15,977 controls (7704 males, 8273 females) of Latino, non-Latino white, or African-American ethnicity remained for analyses. This included a total of 3506 CCRLP controls and 12,471 GERA controls. Case-control association tests were stratified by ethnicity and included 10533 Latino subjects (1949 cases, 8584 controls), 4735 non-Latino white subjects (1184 cases, 3551 controls), and 3972 African-American subjects (130 cases, 3842 controls). SNP associations were calculated using logistic regression, adjusted for the first 10 ancestry-informative principal components generated by Eigenstrat. Single SNP association results were computed assuming an allelic additive model for 0, 1, or 2 copies of the minor allele. Genomic inflation factors were calculated for each ethnicity-stratified GWAS and suggested no substantial inflation of test statistics (λ Latinos = 1.034; λ non-Latino whites = 1.004; λ African Americans = 1.000) following adjustment for the first 10 principal components. The ethnicity-stratified GWAS were subsequently combined using fixedeffects meta-analysis, performed using Plink. SNPs achieving a P value <5.0 × 10 −8 in the meta-analysis were considered to have achieved genome-wide statistical significance by current convention and were carried forward for replication analyses 39 .
Replication datasets. For replication, we used two additional sets of cases and controls. The first is comprised of European-ancestry ALL patients from the COG protocols 9904 and 9905 and controls from the Wellcome Trust Case-Control Consortium. Constitutive DNA from COG samples was extracted from remission blood 40 . DNA samples were genotyped on the Affymetrix 6.0 array and genotype data were downloaded from dbGaP accession phs000638.v1.p1. Genotype data for European-ancestry control samples genotyped on the Affymetrix 6.0 array were downloaded from the Wellcome Trust Case-Control Consortium 41 . Genotyping quality control procedures were conducted independently in cases and in controls 9 . Briefly, samples with genotyping call rates <98% in cases or controls were excluded. SNPs with genotyping call rates <98% were excluded. We excluded subjects showing evidence of non-European ancestry, samples with mismatched reported vs. genotyped sex, and cryptically related subjects (IBD proportion >0.20), and duplicated subjects already included in our discovery set (N = 21). COG and WTCCC genotype data were merged to create a final set of 959 European-ancestry ALL cases and 2624 controls.
A second replication set consisted of 530 Latino ALL cases and 511 Latino controls from the CCLS that did not overlap with either the CCRLP discovery dataset or the COG replication set. Pre-diagnostic constitutive DNA for CCLS samples was extracted from neonatal bloodcards and genotyped using the Illumina Human OmniExpress V1 platform. DNA extraction was performed using the QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA). Samples with genotyping call rates <95% were excluded. Samples were screened for cryptic relatedness using unlinked SNPs and excluded if IBD proportion >0.20. Samples with discordant sex information (reported vs. genotyped sex) were excluded. SNPs with genotyping call rates <95% were excluded. Any SNP with a Hardy-Weinberg equilibrium (HWE) P value <1 × 10 −4 in controls was excluded.
The COG ALL GWAS dataset were obtained from dbGaP study accession phs000638.v1.p1 (Genome-Wide Association Study of Relapse of Childhood Acute Lymphoblastic Leukemia). The ALL Relapse GWAS dataset was generated at St. Jude Children's Research Hospital and by the COG, supported by NIH grants CA142665, CA21765, CA158568, CA156449, CA36401, CA98543, CA114766, CA140729, and U01GM92666, Jeffrey Pride Foundation, the National Childhood Cancer Foundation, and by ALSAC.
Replication association analyses and meta-analysis. Imputation was performed for 1 Mb regions centered on three novel lead SNPs from the discovery meta-analysis. Imputation was performed using the Impute2 v2.1.2 software and its standard Markov chain Monte Carlo algorithm and default settings for targeted imputation 42 . All 1000 Genomes Phase 3 haplotypes were provided as the imputation reference panel 43 . SNPs with imputation quality (info) scores <0.60 or posterior probabilities <0.90 were excluded to remove poorly imputed SNPs. Association statistics for imputed and directly genotyped SNPs were calculated using logistic regression in SNPTESTv2, using an allelic additive model and probabilistic genotype dosages 44 . The effect of individual SNPs on ALL risk was calculated while adjusting for principal components from Eigenstrat.
Association statistics for the COG and CCLS replication datasets were combined using fixed-effects meta-analysis in the META software package to generate a summary odds ratio and P value for the combined replication dataset 45 . The three novel lead SNP associations from discovery analyses were considered to have been successfully replicated if they were statistically significantly associated with ALL risk in this meta-analysis of replication datasets. To preserve the type I error rate, replication P values rejection thresholds were adjusted for multiple comparisons using the Bonferroni-Holm procedure.
A summary meta-analysis of all ALL cases and controls was calculated by combining results from the three ethnicity-stratified discovery samples with the COG and CCLS replication datasets using fixed-effects meta-analysis using the META software 45 . As a measure for between-study heterogeneity, I 2 was calculated 46 .
Allelic functional characterization. Top hit SNPs were associated with genomewide gene expression using two datasets: the Genotype-Tissue Expression (GTEx)