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
] and in vitro
], 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 and , and Figure ). 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 ), 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 and ). 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.
Number of differentially expressed genes by False Discovery Rate (FDR), compared to time 0 hours.
Number of differentially expressed genes by Fold Change (FC), compared to time 0 hours.
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) < (more ...)
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 ). 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%.
Genes regulated 8 hours following dexamethasone treatment.
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 modulate glucocorticoid sensitivity in CEM T-ALL cells [31
]. Other known glucocorticoid response genes upregulated included TSC22D3
] and SOCS1
], both downstream targets of the glucocorticoid receptor, FKBP5
], a co-chaperone of the glucocorticoid receptor, and MAP kinases 5, 6 and 8 [34
]. Downregulated genes at 8 hours included BCL-2
] and C-MYC
], 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 ), and at 48 hours 1 gene was upregulated (t-statistic > 6) and 15 genes downregulated (t-statistic < -6, table ). 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.
Genes regulated 24 hours following dexamethasone treatment.
Genes regulated 48 hours following dexamethasone treatment.
Functional analysis using GSEA and meta-GSEA on the expression profiles obtained at 8 hours and 24 hours after dexamethasone treatment (additional files 1
), 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 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 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 ), which may represent different modes of response to glucocorticoids in patients. Relative to this inter-patient heterogeneity, both cell lines (purple bars, Figure ) and xenografts (black bars, Figure ) 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.
Figure 2 Parametric GSEA of combined top 100 glucocorticoid-induced gene sets with greatest variance from xenograft, patient and cell line models. Hierarchical clustering with gene sets in rows, samples in columns (xenografts - black, patient - green, cell line (more ...)
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 , we proceeded to determine whether we could use fewer replicates and still identify a majority of the differentially expressed genes. Replicate analysis (Figure ) 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.
Figure 3 Recovery scores at 8 hours and 24 hours when randomly selecting all combinations of 3 replicates (3rep) or 2 replicates (2rep) from the set of 4 biological replicates. The Recovery Score represents the proportion of differentially expressed genes from (more ...)