Evaluation of the NOD/SCID xenograft model for glucocorticoid-regulated gene expression in childhood B-cell precursor acute lymphoblastic leukemia

Background Glucocorticoids such as prednisolone and dexamethasone are critical drugs used in multi-agent chemotherapy protocols used to treat acute lymphoblastic leukemia (ALL), and response to glucocorticoids is highly predictive of outcome. The NOD/SCID xenograft mouse model of ALL is a clinically relevant model in which the mice develop a systemic leukemia which retains the fundamental biological characteristics of the original disease. Here we report a study evaluating the NOD/SCID xenograft mouse model to investigate glucocorticoid-induced gene expression. Cells from a glucocorticoid-sensitive xenograft derived from a child with B-cell precursor ALL were inoculated into NOD/SCID mice. When highly engrafted the mice were randomized into groups of 4 to receive dexamethasone 15 mg/kg by intraperitoneal injection or vehicle control. Leukemia cells were harvested from mice spleens at 0, 8, 24 or 48 hours thereafter, and gene expression analyzed on Illumina WG-6_V3 chips, comparing all groups to time 0 hours. Results The 8 hour dexamethasone-treated timepoint had the highest number of significantly differentially expressed genes, with fewer observed at the 24 and 48 hour timepoints, and with minimal changes seen across the time-matched controls. When compared to publicly available datasets of glucocorticoid-induced gene expression from an in vitro cell line study and from an in vivo study of patients with ALL, at the level of pathways, expression changes in the 8 hour xenograft samples showed a similar response to patients treated with glucocorticoids. Replicate analysis revealed that at the 8 hour timepoint, a dataset with high signal and differential expression, using data from 3 replicates instead of 4 resulted in excellent recovery scores of > 0.9. However at other timepoints with less signal very poor recovery scores were obtained with 3 replicates. Conclusions The NOD/SCID xenograft mouse model provides a reproducible experimental system in which to investigate clinically-relevant mechanisms of drug-induced gene regulation in ALL; the 8 hour timepoint provides the highest number of significantly differentially expressed genes; time-matched controls are redundant and excellent recovery scores can be obtained with 3 replicates.


Background
Glucocorticoids such as prednisolone and dexamethasone are critical components of multi-agent chemotherapy protocols used in the treatment of acute lymphoblastic leukemia (ALL) [1] due to their ability to induce apoptosis in lymphoid cells. Despite their use for over 50 years their mechanism of action is not completely understood. Glucocorticoids are steroid hormones that act on target cells through interaction with a specific glucocorticoid receptor (GR) [2]. The GR is held in a cytosolic complex by a number of co-chaperone molecules including HSP-90 and HSP-70 [3], and on ligand binding dissociates from the co-chaperone complex, dimerizes and is transported to the nucleus where it binds to palindromic DNA sequences known as glucocorticoid response elements (GREs) found in the promoter regions of target genes [4]. This leads to the activation of transcription of primary target genes, repression of transcription through interaction with negative GREs [5] or of gene activation through transcription factors such as AP-1 and NF-ΚB [6]. In lymphoid cells, this results in repression of cell cycle progression through cyclin-D3 and C-MYC [7], and cell death through the activation of apoptosis. Glucocorticoids also induce other non-apoptotic mechanisms of programmed cell death including autophagy [8] and mediate a number of pathways involved in the metabolism of carbohydrates, lipids and proteins.
A number of studies have published microarray data of glucocorticoid-induced genes in lymphoid cells, but comparison of the data is complicated by technical differences in platform and chip type. Previous studies of glucocorticoid-induced genes in ALL have been carried out using in vitro cell-line models [9][10][11][12][13][14][15] and patient-derived cells, both in vivo [16] and in vitro [17]. Cell lines are extensively used in the study of ALL but in the process of immortalization acquire multiple genetic defects, particularly in the p53 pathway [18], and mechanisms demonstrated in cell lines are often not replicated in more clinically relevant models. Primary patient cells have a finite supply and rarely survive ex vivo for more than a few days. The nonobese diabetic/severe combined immunodeficient (NOD/ SCID) xenograft mouse model is widely used to study ALL. In this model, human leukemia cells obtained from diagnostic bone marrow biopsies are inoculated into NOD/SCID mice, and on engraftment establish a systemic leukemia which retains the fundamental biological characteristics of the original disease [19]. It has also been shown that the in vivo responses to chemotherapeutic agents, including dexamethasone, correlates with patient outcome [20], and thus the NOD/SCID xenograft mouse model provides a stable, reproducible and clinically relevant model with which to study ALL. Here we report the first study investigating glucocorticoid-induced gene expression in ALL using the NOD/SCID xenograft model, the optimal experimental design, and a comparison of our microarray data to publicly available datasets of glucocorticoidinduced genes in other experimental models.

NOD/SCID xenograft mouse model
All experimental studies were approved by the Human Research Ethics Committee and the Animal Care and Ethics Committee of the University of New South Wales. ALL-3, a glucocorticoid-sensitive xenograft derived from a 12 year old girl with mixed lineage leukemia (MLL)-rearranged BCP-ALL, was chosen for this study. Although MLL-rearranged ALL is associated with a poor prednisolone response and an inferior outcome [21], this patient is currently a long-term survivor. ALL-3 demonstrates in vitro glucocorticoid sensitivity, with an IC 50 of 9.4 nM on exposure to dexamethasone. In the in vivo NOD/SCID xenograft mouse model, ALL-3 is highly responsive to 4 weeks of treatment with single agent dexamethasone, with rapid clearance of leukemic blasts from the peripheral blood and recurrence of leukemia delayed by 63.4 days compared to vehicle-treated controls [20].
Cells from ALL-3 were inoculated by tail-vein injection into 28 NOD/SCID mice. The mice were bled weekly and the samples stained with fluorescein isothiocyanate (FITC)-conjugated anti-murine CD45 and allophycocyanin (APC)-conjugated anti-human CD45 (BioLegend, San Diego, CA). Following lysis of erythrocytes with FACS lysing solution (BD Biosciences, San Jose, CA), samples were analyzed by multiparametric flow cytometry on a FACSCanto cytometer (BD Biosciences, San Jose, CA). Engraftment was calculated as the proportion of human versus total CD45 + cells.
When high level (> 70%) engraftment was achieved in the peripheral blood, between 8 and 10 weeks posttransplantation, the mice were randomized and split into groups of 4 to receive either dexamethasone 15 mg/kg (Sigma-Aldrich, St Louis, MO) or vehicle control by intraperitoneal injection. Mice were culled by CO 2 asphyxiation at 0 hours (pre-treatment, group 1), 8 hours (groups 2 and 3), 24 hours (groups 4 and 5) or 48 hours (groups 6 and 7) following treatment. The mice in groups 6 and 7 received a second dose of dexamethasone or vehicle control at 24 hours. Two mice succumbed early to thymoma, a well-recognized complication in NOD/SCID mice, resulting in 3 mice in each of groups 6 and 7. Cell suspensions of spleens were prepared and mononuclear cells enriched and purified to > 97% human by density gradient centrifugation using LymphoPrep (Axis-Shield, Norway), and cell viability assessed by trypan blue exclusion. RNA was extracted using the RNeasy mini kit (Qiagen, Hilden, Germany) and the RNA integrity verified (Agilent Bioanalyzer, Santa Clara, CA). The RNA was amplified using the Illumina TotalPrep RNA amplification kit (Ambion, Austin, TX) and hybridized onto Illumina WG-6_V3 chips (Illumina, San Diego, CA). The chips were scanned on the Illumina Bead Array Reader (Illumina, San Diego, CA) and gene expression analyzed. The data have been deposited in NCBI's Gene Expression Omnibus [22] and are accessible through GEO Series accession number GSE30392 http://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE30392.

Gene expression and functional analysis
The sample probe profiles with no normalization or background correction were exported from BeadStudio (version 3.0.14, Illumina, San Diego, CA). The data were pre-processed using variance stabilizing transformation [23] and robust spline normalization in lumi [24] which takes advantage of each probe being represented by > 25 beads. Differential gene expression was determined using limma [25] by comparing all treated groups to time 0 hours, with the positive False Discovery Rate correction for multiple testing [26]. Complete linkage hierarchical clustering using Euclidian distance was used to compare groups to each other. Functional analysis was performed using gene set enrichment analysis (GSEA) version 2.04 [27], comparing the limma moderated t-statistic for each probe in a pre-ranked file, against the c2_all collection of gene sets from the Molecular Signatures Database [27]

Comparison of models
The molecular response to glucocorticoids in xenografts was compared to publicly available microarray data [13,16] using parametric analysis of gene set enrichment [28] implemented in the PGSEA package (version 1.20.1, Furge and Dykema) from the Bioconductor project [29], with some modifications to the algorithm to assess significance of the genes that are in the geneset and represented on the microarray, and to allow more control over control sample specification (available upon request). Expression levels of each gene in each sample were converted to expression ratios relative to patient matched controls before glucocorticoid treatment (Schmidt et al), time-matched controls (Rainer et al), or time 0 hours (xenografts). Within each dataset, these gene-level ratios were summarized into geneset-level Zscores, using PGSEA with genesets from the c2_all collection [27]. The Z-scores from each sample from the 3 studies were combined and then compared by hierarchical clustering of the top 100 gene sets demonstrating the greatest variance across the combined studies.

Replicate analysis
The stability of results when reducing the number of replicates was assessed using the Recovery Score method [30] from the GeneSelector package (version 1.4.0) of the Bioconductor project [29].

Results and Discussion
It has been demonstrated that changes in gene expression can be detected as early as 6 hours after treatment of ALL with glucocorticoids in vivo [16] and in vitro [11], although earlier timepoints show few, if any, significantly differentially expressed genes [17]. In this study the 8 hour dexamethasone-treated timepoint demonstrated the highest number of differentially expressed genes compared to baseline control, with far fewer observed at the 24 and 48 hour dexamethasone-treated timepoints (Tables 1 and 2, and Figure 1). Whilst a similar proportion of up-and down-regulated genes were identified at the 8 hour dexamethasone-treated timepoint (1158 vs 1072 respectively, FDR < 0.05), of those with large fold changes (FC > 2 or FC < 0.5, red dots in Figure 1A), 75% were up regulated (199 vs 65 respectively), consistent with the predominant role of glucocorticoids as transcriptional activators. The large numbers of statistically differentially expressed genes (FDR < 0.05) with small fold changes (0.5 < FC < 2) are indicative of both small measurement error across replicates, and thus the high reproducibility of the xenograft model, and good experimental power resulting from using 4 replicates. There was minimal significant differential gene expression across the time-matched controls (Tables 1 and 2). This demonstrates that in the xenograft mouse model, the 8 hour timepoint provides the greatest information, and that these changes are not sustained over later timepoints. The handling of the mice and intraperitoneal vehicle control injections had minimal effect on gene expression, and thus time-matched controls are redundant.
At the 8 hour timepoint, there were 173 genes upregulated with a t-statistic (the ratio of fold change to standard error) > 10 and 25 genes downregulated with a t-statistic < -10 (corresponding to P < 1.74 × 10 -9 and FDR < 2.95 × 10 -7 , table 3). None of these genes had sustained expression changes at 24 or 48 hours, and although this could potentially reflect the early death of sensitive cells, there was no significant difference in the total number of cells harvested from the spleens at any timepoint compared to the time-matched controls (data not shown), and all harvests had a cell viability of ≥ 96%.
The most significantly differentially expressed gene at the 8 hour dexamethasone-treated timepoint was ZBTB16, a known transcriptional repressor and glucocorticoid response gene, which has been shown to , a co-chaperone of the glucocorticoid receptor, and MAP kinases 5, 6 and 8 [34]. Downregulated genes at 8 hours included BCL-2 [35] and C-MYC [36], both previously described, but also HSP90B1, a glucocorticoid receptor co-chaperone and regulator of apoptosis. The only pro-apoptotic gene consistently upregulated across multiple microarray analyses is the BH3-only BCL-2 family member BIM, and it has been shown that BIM has a critical role in glucocorticoid sensitivity and resistance [37], although in this current study BIM was only induced 1.3 fold. Thus these genes identified are consistent with previous reports of glucocorticoid-induced genes in ALL. Within these experimental systems however there are significant potential differences in glucocorticoid exposure between in vitro and in vivo models -a crucial one is that cells in vitro are continuously exposed to glucocorticoid whereas in in vivo models the glucocorticoid is subject to pharmacokinetic and pharmacodynamic changes which more accurately reflect changes in patients. At the later timepoints, significant differential gene expression was much less marked and predominantly downregulated. At 24 hours 5 genes were upregulated (t-statistic > 6) and 10 genes downregulated (t-statistic < -6, table 4), and at 48 hours 1 gene was upregulated (tstatistic > 6) and 15 genes downregulated (t-statistic < -6, table 5). At 24 hours, upregulated genes included NFKBIA, an inhibitor of NF-ΚB, and TRIM74, which was sustained at 48 hours, the significance of which is uncertain. Downregulated genes were those involved in cell cycle progression, including CCNF at 24 hours, and CCNF, CDC20 and AURKA at 48 hours, consistent with growth arrest.
Functional analysis using GSEA and meta-GSEA on the expression profiles obtained at 8 hours and 24 hours after dexamethasone treatment (additional files 1 and 2), revealed a significant upregulation of metabolic pathways, particularly adipogenesis at 8 hours, and a marked effect on pathways associated with cell cycling and proliferation, particularly downregulation of C-MYC at 8 hours and NF-ΚB at 24 hours, and upregulation of Figure 1 Volcano plots of significantly differentially expressed genes following treatment with dexamethasone at 8 hours (A), 24 hours (B), 48 hours (C). Significance was defined as log2 Fold Change > 1 or < -1 with False Discovery Rate (FDR) < 0.05. Each dot represents a single gene, and significant genes indicated by red dots.    apoptotic pathways at 24 hours. Glucocorticoids are known to have effects on multiple cellular metabolic pathways, including glucose and carbohydrate metabolism, and have pro-adipogenic effects [38]. Suppression of C-MYC is a critical step prior to the initiation of apoptosis by dexamethasone in ALL [39] and suppression of NF-ΚB has been described [40].
To determine whether the molecular response to glucocorticoids in this xenograft model of ALL mimicked the effects seen in either glucocorticoid-treated patients with ALL [16] or cell-line models of ALL [13], we applied parametric gene set enrichment analysis (PGSEA) [28]. Comparing gene expression profiles from multiple experiments is notoriously difficult and typically any true similarities are swamped by technical differences in microarray vendor, normalization strategies and analytical approach. By summarizing genes at the gene set level (such as genes in the same pathway), these technical differences are mitigated, allowing comparison of samples from multiple studies.
We performed PGSEA on the 6-8 hour samples from the 3 studies, and then two-dimensional hierarchical clustering to identify the relationships between the different ALL models (Figure 2 and annotated in additional file 3). This revealed considerable heterogeneity in the molecular response to glucocorticoids in patients into at least 2, and possibly 4 different groups (green bars, Figure  2), which may represent different modes of response to glucocorticoids in patients. Relative to this inter-patient heterogeneity, both cell lines (purple bars, Figure 2) and   Figure 2) are remarkably reproducible; we anticipate that adding additional xenograft models of ALL from distinct patients will mirror the heterogeneity of the patient from whom they were derived. It is also evident that overall, glucocorticoid-treated xenografts co-cluster with a group of 3 patients (B-ALL-37, -38, and -40), all of whom had BCP-ALL and a good early prednisolone response, with varied cytogenetics (hyperploidy, t(12;21), and normal respectively). At more relaxed clustering thresholds, the glucocorticoid-treated xenografts cluster with 4 other patients with BCP-ALL (B-ALL-24, -31, -33, and 43) and the cell lines.
We identified 5 clusters of gene sets with distinct expression profiles, each behaving differently in the 3 models of ALL. Cluster 1 demonstrated the markedly heterogeneous patterns seen in patient samples, with the xenograft samples showing a pattern similar to 8 of the patients; cluster 2 showed genesets that showed strong enrichment in the cell line study, and included a number of genesets associated with cell proliferation; cluster 3 did not show any striking differences across the three ALL models; cluster 4 showed genesets downregulated in both xenografts and cell lines compared to the patient samples, and included a number genesets associated with cell cycle progression, DNA/RNA replication and MYC; cluster 5 showed genesets strongly downregulated in the xenograft and cell line models, and included genesets associated with MYC and metabolic processes, particularly catabolism and energy production. In this limited comparison, it is clear that glucocorticoid-induced gene expression patterns seen in ALL are dependent on the experimental model, and that the patterns derived from the xenograft model show a greater similarity to patient-derived data than to cell lines.
A search of the TRANSFAC database v8.3 [41] of CoMoDis [42] identified GRE motifs (within 100 kb either side of the gene) in only 25 (14.5%) of the top 173 upregulated genes at the 8 hour timepoint in this study, and no GRE motifs were identified in the upregulated genes at 24 or 48 hours. This supports accumulating evidence that glucocorticoids exert long-range effects through very distal steroid receptor binding sites [43]. Analysis of significantly differentially expressed glucocorticoid-induced genes in an in vitro cell line study [13] revealed a similar number of early response genes after 6 hours of exposure (60 upregulated (t-stat > 10) and 27 downregulated (t-stat < -10)) but a significantly greater number of genes after 24 hours (593 upregulated (t-stat > 10) and 782 downregulated (t-stat < -10)). Interestingly, all but 2 of the genes upregulated at 6 hours remained significantly upregulated at 24 hours, and 17 of the downregulated genes at 6 hours remained downregulated at 24 hours. GRE motifs were identified in 15 (25.0%) of the top 60 upregulated genes at 6 hours, and 87 (14.6%) of the top 593 genes at 24 hours. The observed difference between the studies in gene expression at later timepoints is consistent with continuous rather than physiological glucocorticoid exposure. In addition, in the cell line study, the GR (NR3C1) undergoes highly significant early and sustained autoupregulation, which in the continuous presence of ligand drives downstream gene expression. In contrast, in the xenograft model minimal GR upregulation is seen at the early timepoint but no significant change in GR expression is seen at either of the later timepoints.
Given the good statistical power observed in Figure  1A, we proceeded to determine whether we could use fewer replicates and still identify a majority of the differentially expressed genes. Replicate analysis ( Figure 3) revealed that at the 8 hour dexamethasone-treated timepoint, a dataset with high signal and differential expression, using data from any 3 randomly chosen biological replicates instead of 4 resulted in excellent recovery scores of > 0.9. That is, on average, 90% of the differentially expressed genes identified from all 4 samples were also identified in any combination of 3 arrays. At 24 hours, a timepoint with less signal, the average recovery score was 0.85 with 3 replicates, but was more variable than at 8 hours. Using just 2 biological replicates recovered 88% of the list of differentially expressed genes at 8 hours, which dropped to 14% at 24 hours. This confirms that the 8 hour time point has the strongest signal, which is reproducible across different subsets of biological replicates. We recommend using a minimum of 3 biological replicates, since fewer replicates destabilized our ability to identify differentially expressed genes. This has important considerations for experimental design, and has significant implications on cost and animal numbers.