Loci on 7p12.2, 10q21.2 and 14q11.2 are associated with risk of childhood acute lymphoblastic leukemia

Acute leukemia is the major pediatric cancer in developed countries, where it affects 30-45 per 10⁶ children per year¹. The childhood acute leukemias are heterogeneous with respect to their underlying cellular and molecular biology, acquired genetic abnormalities and associated clinical responses to combination chemotherapy². It is suspected therefore that the acute myeloblastic leukemia and subtypes of B or T cell precursor acute lymphoblastic leukemia (ALL) may not share a common etiology³.

While epidemiology data are compatible with transplacental carcinogen exposure as a basis for infant leukemias associated with MLL gene fusion³ and dysregulated immune response to common infection is a candidate for childhood ALL³, the role of environmental carcinogenesis in childhood leukemia is presently undefined. It is however probable that the risk of ALL from environmental exposure is influenced by genetic variation. Data from the Swedish family-cancer database lends support to a small familial risk of ALL^4,5, independent of the high concordance in monozygotic twins (which has a non-genetic, in-utero explanation). Although rare (<5% of ALL) direct evidence for an inherited genetic predisposition to ALL is provided by the high risk associated with Bloom’s syndrome, neurofibromatosis, ataxia telangiectasia and constitutional trisomy 21⁶. The heritable basis of susceptibility to ALL outside these syndromes is presently undefined but it is likely that the co-inheritance of multiple low-risk variants contribute to disease risk. Predicated on this hypothesis we have conducted a genome-wide association (GWA) study of ALL analyzing 2 case series genotyped using Illumina Infinium HD Human370 Duo BeadChips.

The first GWA study (GWA-1) was based on genotyping 577 ALL cases from the United Kingdom Childhood Cancer Study (UKCCS)⁷. After applying quality control criteria and exclusion of individuals with non-Western European genotype, SNP genotype data were available for 503 of the cases (Supplementary Fig. 1). For controls we made use of publicly accessible Illumina Hap550K BeadChip genotype data generated on 1,438 individuals from the British 1958 Birth Cohort (58C, also known as the National Child Development Study), which included all live births in England, Wales and Scotland during a single week in 1958⁸.

The second GWA study (GWA-2) was based on genotyping 392 cases from the Medical Research Council (MRC) ALL 97 trial and 36 cases collected by the Northern Institute for Cancer Research (NICR). After applying quality control criteria and exclusion of individuals with non-Western European genotype, SNP genotype data were available for 404 of the cases (Supplementary Fig. 1). For controls we made use of Illumina Hap550K BeadChip genotype data generated on 960 healthy individuals from the UK as part of a study of colorectal cancer (CRC) ⁹.

Across both case series a total of 342,665 tagging SNPs were satisfactorily genotyped (99.7%), with mean individual sample call rates (the percentage of samples for which a genotype was obtained for each SNP) of 99.8%. We excluded 24 individuals because of non-Western European ancestry and 10 because of cryptic relatedness (Supplementary Fig. 2). 293,371 SNP genotypes were available on all 907 cases (824 B-cell ALL, 83 T-cell ALL) and 2,398 controls in the combined data (Supplementary Fig. 1). Prior to undertaking a meta-analysis we searched for potential errors and biases in the 2 GWA studies imposing a high stringency for quality control for SNPs. Specifically, we considered only the 291,473 autosomal SNPs which had call rates >95% in all case and control series, which showed no extreme departure from Hardy-Weinberg equilibrium (HWE; P>10^-5 in controls) and had minor allele frequencies (MAF) exceeding 1% in cases and controls (Supplementary Fig. 1).

Comparison of the observed and expected distributions showed little evidence for an inflation of the test statistics in the 2 datasets (inflation factor¹⁰ λ = 1.034 and 1.002 for GWA-1 and GWA-2 respectively based on the 90% least significant SNPs; Supplementary Fig. 3), thereby excluding the possibility of significant hidden population substructure or differential genotype calling between cases and controls. Using data from both GWA studies we derived joint odds ratios (ORs) and confidence intervals (CI) under a fixed effects model for each SNP and associated P values from the standard normal distribution.

Ten SNPs mapping to 3 genomic regions showed evidence of an association at conventional levels for genome-wide significance (i.e. P < 5 x 10^-7)¹¹. These associations were significant at this threshold irrespective of which control group was used for reference (Supplementary Table 1).

The strongest association signal was attained at 7p12.2 with rs4132601 (combined OR = 1.69, 95% CI: 1.58 – 1.81; P = 1.20 x 10^-19; P_het = 0.68, I² = 0%; Table 1, Supplementary Table 1), which maps to the 3’ region of the ikaros family zinc finger 1 (IKZF1) gene (50,438,098 bps; Fig. 1). The association signal was also highly significant when analysis was confined to B-cell ALL (combined OR = 1.73, 95% CI: 1.61 – 1.85; P = 9.31 x 10^-20; Table 1, Supplementary Table 1). rs6944602 and rs6964823 which map to the 3’ UTR and intron 7 of IKZF1 (50,441,245 bps and 50,427,590 bps; Fig. 1; Supplementary Fig. 4), also displayed statistical support for an association at 7p12.2 at genome-wide threshold (P = 6.02 x 10^-14 and 3.43 x 10^-15 respectively). rs4132601, rs6964823 and rs6944602 map to a 27.3kb block of linkage disequilibrium (LD). LD metrics for rs4132601/rs6964823, rs4132601/rs6944602, rs6964823/rs6944602 are D’ = 1.0, 1.0, 1.0 and r² = 0.42, 0.79, 0.33 respectively. A significant association was also seen with rs7809758, which maps 110kb centromeric to rs4132601 (50,540,827 bps; P = 2.41 x 10^-10), which annotates the Dopa decarboxylase aromatic L-amino acid (DCC) gene. LD between rs7809758 and rs4132601 is not strong (D’ = 0.72, r² = 0.32) and the 9-order difference in statistical support for the association with ALL defined by rs7809758 compared with rs4132601 argues against an independent disease-locus. However, to confirm this will require further mapping and dissection of recombination breakpoints.

Although rs4132601 may not be directly functional the established role of IKZF1 in the biology of ALL strongly implicates variation in IKZF1 as the causal basis of the 7p12.2 association. Ikaros proteins are master regulators of lymphocyte development and differentiation plays a pivotal role in CD4 versus CD8 T-cell lineage commitment decisions¹². Germline mutant mice expressing only non-DNA-binding dominant-negative leukemogenic Ikaros isoforms develop an aggressive form of lymphoblastic leukemia^13,14. Chromosomal deletions involving IKZF1 are common (30%) in high-risk/poor prognosis B-cell precursor ALL¹⁵ and are highly prevalent (95%) in ALL with BCR-ABL1 fusions¹⁶.

To explore the possibility that the association might be mediated through differential IKZF1 expression we investigated the relationship between rs4132601 genotype and mRNA expression level in EBV transformed lymphocytes. Expression was significantly associated with genotype in a dose-dependent fashion (P = 0.005; Fig. 2, Supplementary Fig. 5) with lower levels being associated with risk alleles. This observation is consistent with a model in which the causal variant influences risk by reducing the efficiency of early B-cell differentiation.

The second strongest association signal was attained at 10q21.2 with rs7089424 (combined OR = 1.65, 95% CI: 1.54 – 1.76; P = 6.69 x 10^-19; P_het = 0.29, I² = 7%; Table 1, Supplementary Table 1), which maps to intron 3 of the AT rich interactive domain 5B (ARIDB5) gene (63,422,165 bps; Fig. 1; Supplementary Fig. 4). The association signal was also highly significant when analysis was confined to B-cell ALL (combined OR = 1.70, 95% CI: 1.58 – 1.81; P = 1.41 x 10^-19; Table 1, Supplementary Table 1). Two additional SNPs rs7073837 and rs10740055 having amongst the most extreme P values (P = 4.66 x 10^-16 and 5.35 x 10^-13 respectively) also annotate ARIDB5 (63,369,901bps and 63,388,485bps), localizing to introns 2 and 3 of the gene respectively. These are in LD with rs7089424 (D’ = 0.89, 1.0 and r² = 0.60, 0.43 respectively) and map to a 79.3kb block of LD within ARIDB5. ARIDB5 is member of the AT-rich interaction domain family of transcription factors¹⁷, which plays an important role in embryogenesis and growth regulation¹⁸. While ARIDB5 expression is reported to be upregulated in acute promyelocytic leukemia¹⁹, currently there is no evidence for involvement of ARIDB5 in childhood ALL. Evidence for ARIDB5 having a role in defining B-cell lineage is supported by data from homozygous knockout mice which along with growth retardation phenotype have decreased bone marrow cellularity and reduced numbers of B-cell progenitors¹⁸.

The third strongest statistical evidence for an association was attained at 14q11.2 with rs2239633 (22,658,897 bps; combined OR = 1.34, 95% CI: 1.22 – 1.45; P = 2.88 x 10^-7; P_het = 0.22, I² = 33%; Table 1, Supplementary Table 1). As with the previous regions the association signal was also significant when analysis was confined to B-cell ALL (combined OR = 1.37, 95% CI: 1.26 – 1.49; P = 5.60 x 10^-8; Table 1, Supplementary Table 1). rs2239633 maps within a 25.7kb region of LD which encompasses the gene encoding CCAAT/enhancer-binding protein, epsilon (CEBPE; Fig. 1; Supplementary Fig. 4). Two other SNPs associated with ALL risk at P < 10^-5 (rs7157021 and rs10143875) map within this region of LD, providing additional support for 14q11.2 as a susceptibility locus. CEBP is a suppressor of myeloid leukemogenesis and is mutated in a subset of cases. Intriguingly, CEBPE, along with other CEBP family members has been shown to be occasionally (~1%) targeted by recurrent IGH translocations in B-cell precursor ALL²⁰ suggesting opposing functions of CEBP dysregulation in myeloid and lymphoid leukemogenesis and a possible role in susceptibility to ALL.

Given the biological heterogeneity of ALL, we analyzed the association between the major subtypes of ALL and rs4132601, rs7089424 and rs2239633 genotypes through case-only logistic regression (Table 2). The primary impact of variation defined by the 7p12.2, 10q21.2 and 14q11.2 risk variants is for B-lineage leukemia. Subtype analysis of B precursor ALL provides strong evidence that variation at 10q21.2-ARIDB5 is highly associated with the risk of developing hyperdiploid ALL (P = 3.84 x 10^-6; Table 2).

Fine-mapping and resequencing is required to identify the specific functional variant underlying each of the associations we have identified. Accepting HapMap is not comprehensive; few non-synonymous SNPs have been documented in IKZF1, ARID5B, and CEBPE and none are correlated with rs4132601, rs7089424 or rs2239633. These data suggest the associations identified are mediated through LD with sequence changes that influence gene expression rather than protein sequence or through LD with low frequency variants that are not catalogued by HapMap. While we did not find a significant relationship between genotypes and ARID5B and CEBPE expression (Supplementary Fig. 5) we were able to demonstrate a relationship between rs4132601 genotype and IKZF1 expression compatible with the causal variant at this locus influencing differential expression.

When we modeled pairwise combinations of the SNPs, we did not find evidence of interactive effects between any of the three loci identified (P > 0.1 for all pairwise interactions: Supplementary Table 2), suggesting that each locus has an independent role in ALL development. While the risks of ALL associated with 7p12.2, 10q21.2 and 14q11.2 variants are modest the carrier frequencies of risk alleles of rs4132601, rs7089424 and rs2239633 are high in the European population and hence the loci make a significant contribution to the development of ALL, underlying ~64% of cases.

Our findings provide the first unambiguous evidence that common genetic variation influences the risk of developing pediatric ALL and a strong rationale for searching for additional risk variants through additional GWA scans. Furthermore, these findings provide novel insight into the development of ALL. It is striking that the 3 risk variants we identify map to genes involved in transcriptional regulation and differentiation of B-cell progenitors. Ethnic differences in the risk of ALL are well recognized¹, it will therefore be interesting to explore how our findings translate to non-Western European populations.