Integrated genomics of susceptibility to alkylator-induced leukemia in mice

Therapy-related acute myeloid leukemia (t-AML) is a secondary malignancy attributable to chemotherapy and/or radiation exposure. t-AML comprises 5-20% of adult AML cases and its prevalence is increasing along with the size of the population undergoing chemotherapy [1, 2]. While chemotherapy regimen [3] and genetic background [4] contribute to t-AML, the risk factors are not well understood. Strong evidence for genetic predisposition to t-AML is provided by inherited cancer syndromes such as neurofibromatosis, where germline mutations of NF1 are associated with increased risk of t-AML in humans and mice [5, 6]. Gaining a better understanding of t-AML susceptibility factors is a pressing concern as it may lead to prevention strategies and provide insight into the genesis of de novo AML.

One class of chemotherapeutics associated with t-AML is the alkylators (i.e. melphalan, busulfan, thiotepa). The therapeutic effect of alkylator agents is believed to result from the formation of DNA adducts and single and double-strand breaks, which trigger apoptosis or growth arrest [7]. Based on this presumed mechanism of alkylator action, genes involved in DNA repair [8], response to oxidative stress [9], and drug metabolism [10] have been investigated as mediators of t-AML susceptibility in candidate gene studies, with largely inconclusive results. A recent study in our lab investigated the genetic basis of t-AML susceptibility using inbred mice [11]. In this study, eight to twelve individual mice from each of 20 inbred strains were treated with the alkylating agent N-nitroso-N-ethylurea (ENU), a potent mutagen with a propensity to cause AT:TA transversions and AT:GC transitions [12]. Mice were monitored for the development of AML for up to 16 months post ENU exposure. The incidence of AML varied by strain from 0 to 80% (H² = 0.10, P-value < 0.001), supporting the hypothesis that there is a strong genetic component in t-AML susceptibility.

We hypothesized that the pre-exposure transcriptional state of hematopoietic stem and progenitor cells, the putative target of leukemogenesis [13], underlies variation in susceptibility to t-AML. A pre-exposure transcriptional basis of susceptibility would be expected if a rapid response is critical in determining a cell's ultimate fate upon mutagen exposure. This hypothesis is consistent with the observation that expression of genes critical to surviving genotoxic stress in yeast does not change after exposure to DNA-damaging agents [14], implying that the necessary factors are already expressed at baseline. A similar situation has been reported in human lymphoblastoid cell lines, in which the pre-exposure transcriptional state of the cell more accurately predicts survival after alkylator treatment than the post-exposure state [15].

In this study, we apply an integrated genomics approach [16] to identify and prioritize genetic and transcriptional networks underlying t-AML susceptibility in mice (Figure 1). By linking expression profiles and complex traits to common genomic loci, this method can ameliorate some of the limitations inherent in genetic association and expression profiling studies [17–21]. When combined with network analysis, this methodology has proven useful in elucidating the biological pathways underlying several complex traits [22, 23].

There is accumulating evidence that many genetic contributors to complex traits are not protein-coding changes [32, 33]. If true, then other classes of genetic events that can affect phenotype must, at some level, impact gene expression (i.e., eQTLs). Hypothesizing that such events contribute to t-AML susceptibility, we took an integrated genomics approach to identify and prioritize candidate transcriptional networks. The first step in this approach was to identify eQTLs in hematopoietic stem and progenitor cells, the likely target of leukemic transformation. Previously, we described a CNVR eQTL map in classical inbred mice [24]. In the current work, we expanded this map to include SNP-based haplotype eQTLs. In deriving the mouse haplotype map, we found surprisingly little correlation between haplotypes and neighboring CNVRs. This is in contrast to human studies, where nearly 75% of common CNVRs are estimated to be in linkage disequilibrium with neighboring SNPs [34]. This suggests that at the currently available resolution and coverage (and genotyping accuracy), mouse haplotypes and CNVRs represent distinct sources of genetic information. We found two-fold more CNVR eQTLs than haplotype-based eQTLs (401 vs. 167). It is tempting to speculate that this difference in eQTL types is because CNVRs have a stronger impact on expression in cis and therefore are more likely to be detected as eQTLs. However, the difference could largely be due to the reduced power to detect haplotype eQTLs because of the exacerbated multiple testing problem that comes with performing approximately 20 times more statistical tests. In total, greater than 60% of the eQTLs were reproducible in an independent dataset.

The second step in the integrated approach was to find genes differentially expressed between t-AML susceptible and resistant strains. Because unsupervised clustering using all expressed transcripts grouped strains by susceptibility status, we expected to find a large number of genes associated with susceptibility. Greater than 7% of the expressed transcripts are differentially expressed (976/13,496). These genes are enriched in several, independent biological processes, most notably apoptosis. Among the differentially expressed intrinsic apoptosis genes are Caspase 9 (Casp9), B-cell leukemia/lymphoma 2 (Bcl-2), BCL2-associated agonist of cell death (Bad), BCL2-associated X protein (Bax), and mutS homolog 6 (Msh6). Msh6 is a member of the mutSα DNA mismatch recognition complex that has been shown to mediate apoptosis in certain contexts [35, 36]. Notably, the absence of mutSα activity in myeloid progenitors results in the complete loss of O 6-methylguanine (O 6MeG)-mediated cytotoxicity [37]. That resistant strains have higher expression of Msh6 suggests that upon alkylator exposure, resistant strains may recognize DNA damage and respond appropriately (i.e. undergo apoptosis) whereas the KL cells of susceptible strains may tend to survive, accumulate mutations, and become transformed. In KL cells, almost all susceptible strains have no detectable expression of Casp9, an initiator of programmed cell death, suggesting that these cells (low-to-no Casp9 expression) are less primed for Casp9-dependent apoptosis. Although knockout of Casp9 in mice is embryonic lethal, presumably due to severe brain development defects [38, 39], Casp9 is not required for apoptosis during normal hematopoiesis [40]. Therefore, inbred strains may exhibit genetically-driven yet tissue-specific differences in apoptosis such that hematopoietic cells of susceptible strains are relatively protected from cell death after exposure to genotoxic stress.

Differential expression and gene enrichment analysis highlighted several biologically plausible pathways that may underlie t-AML susceptibility. However, it remained unclear which pathway members, if any, are causal contributors to the phenotype, as illustrated by the complex expression patterns of the intrinsic apoptosis genes. More broadly, the role and relative importance of each of the 917 differentially expressed genes in susceptibility remained undetermined. We hypothesized that among the 917 differentially expressed genes would be a subset in which expression variation is caused by cis-encoded genetic variation. Further, we posited that these 'anchor' genes cause expression variation of multiple downstream genes, which collectively are associated with t-AML susceptibility (or resistance). The mechanisms by which anchors might act in trans are varied. They include altered transcription factor abundance and homeostatic or compensatory forces within and between biological pathways. Regardless of the mechanisms of action, the identification of putative events that influence susceptibility and their linkage to gene networks forms a powerful and practical strategy to both find biological pathways underlying cancer susceptibility and to prioritize candidate mediators. Therefore, as the third step in the integrated genomics approach, we identified networks of genes that are significantly correlated with candidate susceptibility anchors. To validate the networks, we used independent gene expression data from multiple hematopoietic populations, trimming the networks of response genes whose expression was not reproduced.

One of the benefits of the integrated genomics approach is that it can implicate biological processes that would not have been detected using differential expression alone. The susceptibility networks that we identified are enriched in genes involved in DNA repair, base excision repair, apoptosis, and cell cycle, among other annotations. A second potential benefit of the integrated approach is that it differentiates between upstream (anchors) and response genes, an advantage over existing approaches that derive gene regulatory networks from expression data alone. While this is a hypothesis that remains to be tested, the identification of candidate upstream factors will be useful in prioritizing among apoptosis-related genes for experimental validation. For example, although Casp9 and Bcl2 are differentially expressed, Casp9 is also the candidate anchor of module A_33, the module most strongly associated with susceptibility status. We speculate that perturbation of candidate anchors, such as Casp9, are more likely to be informative in elucidating susceptibility than response genes (i.e. Bcl2).

Network analysis allowed us to predict the function of uncharacterized genes. For example, A630001G21Rik is expressed primarily in primitive hematopoietic and B-cells [41], yet its function is undetermined. Our analysis places it as the anchor of module A_12, which is enriched in apoptosis-related genes, including Bcl2. Therefore, A630001G21Rik may play previously unknown role in regulation of Bcl2 expression and apoptosis activity. Similarly, Cytoskeleton-associated protein-like 2 (Ckap2l) is the anchor of the largest module, A_16, enriched in both cell cycle and DNA repair genes. Although Ckap2l is highly expressed in hematopoietic progenitors [41], its functions are unknown. Its closest ortholog, Ckap2, is highly expressed in mouse stem cell lines and has detectable expression in hematopoietic progenitors, bone marrow, osteoclasts, osteoblasts, and macrophages [41]. There is a growing body of literature suggesting that Ckap2 (also known as Tumor-associated microtubule-associated protein) is involved in cell cycle progression [42–44]. It is possible that Ckap2l contributes to cell cycle regulation in HSCs and progenitors, and that genetic disturbances of its expression impact t-AML susceptibility. Experiments that perturb expression of anchor genes such as Casp9 and Ckap2l to assess their impact on module expression and activity are the next logical steps in determining the role of candidate networks in susceptibility. If such experiments demonstrate a causal link between anchor genes and module expression, then moving forward to formally define their role in leukemia will be warranted.

A drawback to the anchored network approach, as currently implemented, is that it assumes there is only a single anchor per module. In cases where CNVRs disrupt regulatory elements, it is possible that a single genetic event impacts the expression of multiple neighboring genes. For example, in module A_37 (Figure 4C) we found 10 response genes within 7 Mb of a CNVR. This module warrants special attention because it includes poly (ADP-ribose) polymerase family member 2 (Parp2, the anchor) and apurinic/apyrimidinic endonuclease 1 (Apex1), both members of the base excision repair pathway [45, 46]. Both genes have lower expression in susceptible strains, again suggesting that lowered overall DNA damage response promotes susceptibility.

A caveat to the current work is that maps of genetic variation in the mouse genome are incomplete, a knowledge gap that promises to be filled by more informative SNP arrays [47] and next-generation sequencing [48, 49]. It is possible that un-captured genetic variants may be the ultimate cause of the observed co-expression networks. These variants may mediate their impact through mechanisms other than altering the expression of anchors. In the extreme case, all modules may not be controlled by anchor expression, but by undetected causes. Nevertheless, the modules themselves are still informative in that they describe sets of coordinately regulated genes that, collectively, are associated with both t-AML susceptibility and biologically plausible processes and pathways.

To our knowledge, this is the first report of an integrated genomics approach to dissect the role of the pre-exposure transcriptional state in t-AML susceptibility. From a clinical perspective, t-AML is important because the response to treatment is poor and survival is short [3]. But because t-AML is a clinically-induced malignancy, it is by definition preventable. Therefore, a long-term goal in this field is to gain sufficient understanding of susceptibility factors in order to make worthwhile the personalization of chemotherapeutic regimens based on t-AML risk. The transcriptional networks and their candidate anchors described here are an important early step towards gaining such an understanding.