Toxicogenomics profiling of bone marrow from rats treated with topotecan in combination with oxaliplatin: a mechanistic strategy to inform combination toxicity

Combinations of anticancer agents may have synergistic anti-tumor effects, but enhanced hematological toxicity often limit their clinical use. We examined whether “microarray profiles” could be used to compare early molecular responses following a single dose of agents administered individually with that of the agents administered in a combination. We compared the mRNA responses within bone marrow of Sprague-Dawley rats after a single 30 min treatment with topotecan at 4.7 mg/kg or oxaliplatin at 15 mg/kg alone to that of sequentially administered combination therapy or vehicle control for 1, 6, and 24 h. We also examined the histopathology of the bone marrow following all treatments. Drug-related histopathological lesions were limited to bone marrow hypocellularity for animals dosed with either agent alone or in combination. Lesions had an earlier onset and higher incidence for animals given topotecan alone or in combination with oxaliplatin. Severity increased from mild to moderate when topotecan was administered prior to oxaliplatin compared with administering oxaliplatin first. Notably, six patterns of co-expressed genes were detected at the 1 h time point that indicate regulatory expression of genes that are dependent on the order of the administration. These results suggest alterations in histone biology, chromatin remodeling, DNA repair, bone regeneration, and respiratory and oxidative phosphorylation are among the prominent pathways modulated in bone marrow from animals treated with an oxaliplatin/topotecan combination. These data also demonstrate the potential for early mRNA patterns derived from target organs of toxicity to inform toxicological risk and molecular mechanisms for agents given in combination.


INTRODUCTION
Effective anticancer treatments generally require the use of a combination of drugs. A challenge for combination treatment is a shift in the extent and severity of adverse effects for the combination compared with either of the agents when given alone. Most Phase I trials for targeted agent combinations have been designed to escalate the drug doses to the maximally tolerated dose (MTD). Data that are now available suggest that combinations of molecularly targeted agents (MTAs) can lead to more severe toxicities and that dose reduction is often required for regimens to be tolerable (Kummar et al., 2010). Although empirical approaches are currently being used, a preclinical assessment that can provide a mechanistic basis for adverse effects of drug combinations is desirable. Systems-based computational approaches that integrate, transcriptomics, proteomics, or other profiling data may be useful to discover network-level alterations in normal cell signaling in normal tissues after in vivo administration.
Toxicogenomics, the application of transcription profiling to toxicology (Nuwaysir et al., 1999), has been widely used for elucidating the molecular and cellular actions of drugs and chemicals on biological systems, flagging potential for toxicity (before functional damages occur), and providing classification of known or new toxicants based on signatures of gene expression (Guerreiro et al., 2003;Ruepp et al., 2005;Bushel et al., 2007;Chengalvala et al., 2007;Ganter et al., 2008;Ryan et al., 2008). In the current study, we examined the benefit of applying transcriptomics to assess risk of enhanced tar-get organ toxicity when two drugs are used in combination. We hypothesized that an additive toxicity pattern can be inferred from comparative analysis of early mRNA responses of tissues obtained following a single dose of two agents given individually with the agents given in combination. We also anticipate that biological pathways revealed by gene response patterns will provide a more comprehensive understanding of toxicity and provide a mechanistic basis for experimental investigation.
For this initial proof of concept, we used the chemotherapeutic agents topotecan and oxaliplatin, which have been explored for use in combination for treating various cancers (Alexandre et al., 2002;Tortora et al., 2002;Elkas et al., 2007). Topotecan is a topoisomerase I (Topo I) inhibitor that forms a complex on DNA leading to double-stranded DNA breaks and ultimately cell death. Oxaliplatin is a diaminocyclohexane platinum compound that acts as a DNA alkylating agent. Pre-exposure to oxaliplatin transiently increased Topo I-mediated single-strand breaks, suggesting that DNA platination might stimulate Topo I DNA cleavage activity providing a rationale for the use of platinum-based compounds with topoisomerase inhibitors (Goldwasser et al., 1999). Combinations of topo I inhibitors and platinum derivatives, in general, have synergistic anti-tumoral effects, but their clinical use is limited by hematological toxicity, which is somewhat dependent on the sequence of drug administration. We also wanted to examine the possibility that the acute response of a target organ to a single dose of treatment and obtaining an associated molecular footprint, can aid in investigation about mechanisms of target-mediated toxicity. In our study, a single dose of topotecan was given to Sprague Dawley rats either before or after administration of oxaliplatin to determine whether combination administration or the order of administration influenced bone marrow toxicity and the microarray signature profile obtained following combined administration of the drugs. We found that a single dose intravenous infusion of topotecan or oxaliplatin, given alone or in combination, to male rats resulted in the formation of histopathological lesions in bone marrow (hypocellularity) observed after only1 h of topotecan administration or 6 h after oxaliplatin administration. The severity of bone marrow lesions increased when topotecan was given prior to oxaliplatin compared with oxaliplatin given prior to topotecan indicating that severity of toxicity was affected by sequence of administration.
We also characterized the molecular and pathway signatures in bone marrow for a topotecan/oxaliplatin combination based on global gene expression analyses and comprehensive bioinformatics profiling. Based on our results, we propose that single dose rodent studies and microarray analysis of mRNA patterns derived from bone marrow represents a mechanistic approach to evaluate the potential risk of enhanced toxicity for combination agents.

TEST ARTICLES
Oxaliplatin (1.2 g; NSC 266046) and topotecan (310 mg; NSC 609699) were supplied by the National Cancer Institute (NCI) and received at Southern Research Institute on October 19, 2010. The test articles were received on dry ice and then refrigerated. The following reagents were received from the NCI chemical repository and stored at room temperature and used in the preparation of dose formulations: 5% dextrose injection, USP (PSS World Medical; Kennesaw, GA); Saline Solution 0.9% (saline; Nova Tech Inc.; Grand Island, NE).

DOSE FORMULATION PREPARATION
Dose formulations of topotecan were prepared in saline to contain a nominal concentration of topotecan hydrochloride of 1.2 mg/mL. For preparation, the appropriate amount of saline was added to the required amount of topotecan hydrochloride and the formulation was stirred until a solution was obtained. Dose formulations of oxaliplatin were prepared in 5% dextrose in water (D5W) to contain a nominal concentration of oxaliplatin of 3.75 mg/mL. For preparation, the appropriate amount of D5W was added to the required amount of oxaliplatin and the formulation was stirred until a solution was obtained.

ANIMALS AND TREATMENTS
The 96 Sprague-Dawley male rats used in this study were obtained from Charles River Laboratories, Inc. (Raleigh, NC). Each rat was procured with an indwelling femoral vein cannula. Animals were given a unique identification number by ear punch. Prior to dosing on day 1, the catheter of each rat was checked for patency. Rats with patent catheters were randomly assigned to one of 8 treatment groups. On day 1, the rats were approximately 8 weeks of age and weighed between 271.2 and 342.2 g. Teklad Certified Rodent Diet 2016C (Harlan; Madison, WI) and tap water (Birmingham public water supply) were provided ad libitum to the rats prior to and throughout the study. The animals were individually housed in solid-bottom polycarbonate cages on stainless steel racks in a room maintained at a temperature of 70-79 • F and a relative humidity of 44-66%. Room lights were controlled by an automatic timer set to provide a 12/12 light:dark cycle. Heat-treated hardwood chip bedding (P. J. Murphy Forest Products, Corp.; Montville, NJ) was used as bedding material. No known contaminants were present in the food, water, or bedding that would be expected to interfere with or affect the outcome of the study. Cage size and animal care conformed to the guidelines of the Guide for the Care and Use of Laboratory Animals (National Research Council, 2011) and the U.S. Department of Agriculture through the Animal Welfare Act .
The study design (with doses administered) is provided in Table 1. A dose slightly below the MTD was chosen for these studies. Since the goal of our study was to investigate gene responses induced by a single dose that could potentially be premonitory for toxicity, we administered a dose that would not be confounded by marked pathological changes in the tissue. The MTD for oxaliplatin in rats was referenced in the Center for Drug Evaluation and Research pharmacology review of NDA # 21-492 (FDA, 2002) and is close to the human clinical dose of 85 mg/m 2 . The dose selected for topotecan was derived from the time course of the hematological effects of topotecan after two consecutive daily administrations to tumor-bearing rats that was described using a semiphysiological model to predict the hematotoxic effects of topotecan (Friberg et al., 2002;Segura et al., 2004). Each rat was administered a single intravenous (IV) infusion dose, given over a 30 min interval, of topotecan (of 4.7 mg/kg), oxaliplatin (15 mg/kg), or respective vehicle control formulation (see Table 1). Thirty minutes after the end of the first infusion, topotecan, oxaliplatin, or a respective vehicle control formulation was administered by IV infusion, given over a 30 min interval. For example, animals given topotecan (or corresponding vehicle) during the first infusion were given oxaliplatin (or corresponding vehicle) during the second infusion. All doses were delivered at a flow rate of approximately 4 mL/kg/30 min. Infusion volumes were based on the mean body weight of animals on the day prior to dosing. Rats were observed twice daily during the pre-study and study periods for signs of mortality and morbidity. Detailed clinical examinations of each rat were collected prior to euthanasia. Each animal was weighed on the day of dosing prior to dose administration (for the calculation of infusion rate; data not reported).

HISTOPATHOLOGY
Samples of femur bone marrow were collected for histopathology and fixed in 10% neutral buffered formalin. All slides were examined by a veterinary pathologist. Each lesion was listed and coded by the most specific topographic and morphologic diagnoses, severity, and distribution using Toxicology Data Management System (TDMS) nomenclature. A four-step grading system was used to rank the severity of microscopic lesions for comparison among groups.

TISSUE COLLECTION AND RNA ISOLATION
At 1, 6, or 24 h after the completion of dosing, 4 rats per group were euthanized by CO 2 asphyxiation. Immediately after euthanasia, bone marrow from left femurs was collected for RNA isolation. Bone marrow was flushed from bone using RNAlater® and refrigerated (approximately 2-8 • C) for at least 24 h and then stored at or below −20 • C. RNA isolation from individual tissues was accomplished using an RNeasy microarray kit (QIAGEN Inc., Valencia, CA). After extraction, the RNA concentration of each sample was determined using a RiboGreen assay. The RNA Integrity Number (RIN) of each sample was determined using an Agilent 2100 Bioanalyzer with 2100 Expert Software (Version B.02.06.S1418; Agilent, Santa Clara, CA). Only samples with a RIN of 7.5 or higher were deemed acceptable for gene expression analysis. Two samples, one from group 3 and a group 2 did not meet the RIN for further processing. Each of the remaining RNA samples were diluted to the required concentration (250 ng/μL) and then stored at or below −70 • C prior to shipment to Expression Analysis (Durham, NC) for microarray hybridization.

MICROARRAY HYBRIDIZATION
Microarray data was generated by Expression Analysis (Durham, NC). RNA samples were converted into labeled target antisense RNA (cRNA) using the Single-Round RNA Amplification and Biotin Labeling System (Enzo Life Sciences, Farmingdale, NY). Briefly, 2.5 μg of total RNA was converted into double stranded cDNA via reverse transcription using an oligo-d(T) primer-adaptor. This cDNA was purified and used as a template for in vitro transcription using T7 RNA polymerase and biotinylated ribonucleotides. The resulting cRNA was purified using magnetic beads and quantitated using spectrophotometry. Next, 11 μg of purified cRNA was fragmented using a 5X fragmentation buffer (200 mM Tris-Acetate, pH 8.1, 500 mM KOAc, 150 mM MgOA), then a hybridization cocktail was prepared and added to the fragmentation product using the Hybridization, Wash and Stain kit (Affymetrix, Santa Clara, CA), applied to Affymetrix GeneChip Rat Genome 230 2.0 Arrays, and incubated at 45 • C for 16 h. Following hybridization, arrays were washed and stained using standard Affymetrix procedures before scanning on the Affymetrix GeneChip Scanner 3000 using factory PMT settings. Data extraction was completed with Expression Console software using a target scaling of 500. The data are available in the Chemical Effects in Biological Systems (CEBS) repository (Waters et al., 2008) under investigation accession number: 007-00005-0000-000-2 and study accession number: 007-00005-0001-000-3 as well as in the Gene Expression Omnibus (GEO) (Edgar et al., 2002;Barrett et al., 2013) under accession number GSE63902.

DETECTION OF DIFFERENTIALLY EXPRESSED GENES (DEGs)
Data analysis was performed using the Bioconductor R package (Gentleman, 2005). We used the package "affyQCReport," combined with principal component analysis (PCA), for array quality control (i.e., outliers detection). After the one outlier sample was excluded (sample 1344-05H, 6 h oxaliplatin followed by vehicle), the remaining Affymetrix raw CEL files were preprocessed using the robust multichip average (RMA) algorithm (Irizarry et al., 2003a,b); which includes background correction, quantile normalization, and summarization by the median polish approach. Next, we performed PCA-based gene filtering on the log base 2 scale data from RMA using the package "pvac," where the filtering is based on a score measuring consistency among probes within a probe set (Lu et al., 2011). Finally, the statistical significance of DEGs was accessed by comparison of treated samples to timematched controls using the empirical Bayes based method known as a limma t-test, which is available in the bioconductor package limma (Smyth, 2005). Analyses were carried out in two batches, where the first batch includes samples from groups 1, 3, 5, and 7, and the second batch contains samples from groups 2, 4, 6, and 8, (see study design Table 1).

DETECTION OF PATTERNS OF CO-EXPRESSED GENES
log base 2 scaled data of the probes remaining after the PCAbased filtering were analyzed for expression patterns across the samples using methodology known as extracting gene expression patterns and identifying co-expressed genes (EPIG) . For each treatment group, the pixel intensity data for each probe was converted to a ratio value by dividing the average probe pixel intensity by the 1 h control samples from that group and then taking the log base 2. EPIG uses correlations (r) across all the sample groups, signal to noise ratio (s/n) within groups of samples, and magnitude of fold change (FC) for a probe within a group to first detect all potential patterns in the data, and then categorize each probe to the pattern that is most statistically significant in terms of the correlation between the probe profile and the pattern. The parameter settings for the EPIG analysis were the defaults: r = 0.8, s/n = 2.5, and FC = 0.5. We used a minimum pattern cluster size of 6 for finding all potential patterns.

ENRICHMENT OF BIOLOGICAL PROCESSES BASED ON DEGs
For each DEGs list derived at p < 0.01 and absolute fold change >2, the Affymetrix GeneChip Rat Genome 230 2.0 Array probe sets were mapped to the Gene Ontology (GO) biological processes (BPs) of the genes they represent using version 2.10.1 of the GO database and the rat2302 database. The 2.14.0 version of the "topGO" Bioconductor package in R was used to perform enrichment of GO BP terms. The classic algorithm (where the significance of a node is considered independent of the significance of the neighboring nodes) and the Kolmogorov-Smirnov test statistic were used for enrichment with node size of at least 5 DEGs. Significance of enrichment was set to p < 0.05. The union of the enriched GO BPs terms from all of the DEGs lists yielded 641 terms. BP terms which were not significant had missing pvalues and were imputed with 1.0. The 641 GO BPs terms were clustered based on the −log base 10 p-values and using Pearson correlation (r) as the dissimilarity metric with average linkage grouping. Clusters (those with a 1-r ≥ 0.9) of GO BP terms (n ≥ 30) were labeled according to the node having the maximum number of paths to it within the GO BP subtree directed acyclic graph derived from the terms in the cluster.

PATHWAY ANALYSIS OF GENES WITHIN PATTERNS EXTRACTED BY EPIG
Clusters of genes identified by EPIG were analyzed using Pathway Studio 9® (Ariadne Genomics, Rockville, MD), to find enriched pathways. Enrichment analysis in Pathway Studio 9® was performed by Gene Set Enrichment Analysis (GSEA) and Sub-Network Enrichment Analysis (SNEA) algorithms. Functional enrichment was performed using Fisher's exact test. SNEA enrichment in Pathway Studio was calculated using the Mann-Whitney test, a non-parametric method for comparing the medians of two distributions. Significant enrichment was set at p < 0.05. For GSEA, pathways with fewer than 3 entities represented were filtered from the data sets.

VISUALIZATION OF GENE EXPRESSION ON PATHWAYS
To visualize the changes of gene expression within a particular pathway, Pathvisio version 3.1.3 (Van Iersel et al., 2008) was used to overlay the average of the log base 2 ratio values from the replicate samples onto biological pathways obtained from Wikipathways (Kelder et al., 2012). The cluster IDs from UniGene (according to the rat Rn4 July 12, 2012 release of refSeq version 54) was used to map the probe sets on the microarray to the genes on the Wikipathways.

GENE REGULATORY NETWORK RECONSTRUCTION
The Gene Regulatory Network Inference (GRNInfer) software (Wang et al., 2006) with default parameter settings (λ = 0.0 and threshold = 1 × 10 −6 controlling the sparseness and the complexity of the network respectively) was used to reconstruct the interactions of the 19 genes in the p53 signaling Wikipathway that had gene expression data mapped. The gene expression data from microarray probes mapping to the same UniGene cluster were averaged and then biological replicates at each time point were averaged. GRNInfer uses linear programming and singular value decomposition of the time point gene expression data to derive of the interactions of the genes within the network given the above thresholds. For topotecan followed by oxaliplatin and oxaliplatin followed by topotecan, the averaged gene expression measures at each time point were used to reconstruct the gene regulatory network depending of the order of administration of the two drugs.

CLINICAL OBSERVATIONS
No mortality occurred during the study and no drug-or treatment-related adverse clinical signs were observed for any animal during the study. The mean body weight of animals in the individual dose groups ranged from 300 to 310 g on the day of dosing.

HISTOPATHOLOGY
Bone marrow hypocellularity was observed in the bone marrow of animals dosed with topotecan, oxaliplatin, or a combination of these two drugs (Figure 1 and Table 2). Minimal bone marrow hypocellularity consisted of an approximately 10-20% decrease in the normal population of cells which reside in the bone marrow in comparison with control animals. The remaining cell population consisted primarily of band and mature granulocytes (neutrophils and eosinophils), mature erythrocytes, and megakaryocytes. Occasionally apoptotic cells were observed in the bone marrow of some animals. These lesions were observed as early as 1 h after the end of topotecan administration and 6 h after the end of oxaliplatin administration. Bone marrow hypocellularity was observed at a higher and earlier onset of incidence in animals that received topotecan alone or topotecan in combination with oxaliplatin than in dose groups that received only oxaliplatin. One-hour post-treatment, bone marrow hypocellularity was observed for animals in all of the dose groups administered topotecan alone or topotecan in combination with oxaliplatin, but not in the dose groups given only oxaliplatin. There also appeared to be an increase in severity of bone marrow suppression from mild to moderate when topotecan was given prior to oxaliplatin than when oxaliplatin was given prior to topotecan. These data support a sequence dependency for the severity of bone marrow toxicity for this combination.

RNA ISOLATION
The integrity of each isolated RNA sample was equal to or greater than the protocol-specified minimum RIN (7.5), with the exception of two samples. Acceptance criteria could not be obtained for bone marrow collected from one of the animals from group 3 (topotecan/D5W treated group) and one animal in group 2 (D5W/saline vehicle control group). These samples were not submitted for gene expression analysis.

PATTERNS OF CO-EXPRESSED GENES
We focused our analysis on uncovering groups of genes displaying similarity in their expression patterns and comparing the gene responses between treatment groups. The EPIG approach utilizes the underlying structure of gene expression data to extract patterns and identify co-expressed genes that are responsive to experimental conditions . The response patterns of genes were used to select gene sets that represent potential signatures for effects of the combination treatments that differed from that of the single agents. As shown in Figures  Minimal hypocellularity in male rat exposed to oxaliplatin. (C) Mild hypocellularity in male rat exposed to topotecan. (D) Moderate hypocellularity in male rat exposed to oxaliplatin followed by topotecan. Bone marrow hypocellularity was graded based on the estimated percentage of cell loss with minimal hypocellularity = 10-20% cell loss, mild hypocellularity = 30-40% cell loss, moderate hypocellularity = 50-60% cell loss, and marked hypocellularity = greater than 60% cell loss.

www.frontiersin.org
February 2015 | Volume 6 | Article 14 | 5 respectively. Patterns that identified differences between individual treatment and combination treatment groups are marked with an asterisk in Supplementary Figures 1-3. The Venn diagrams in Supplementary Figure 4 reveal that there is little overlap of co-expressed genes changes (of 4-fold or more) between oxaliplatin and topotecan except for at the 24 h time point where 73 genes overlap. The heat map in Figure 3A shows gene expression patterns across all treatment groups and time points. There is a clear time response in the patterns of gene expression from the treatments where by the 24 h time point, the patterns of expression are either maximally induced or maximally repressed. A list of representative pattern-specific genes showing at least 4-fold differences from their respective controls at the 24 h bone marrow collection time point is provided in Table 3 (up-regulated) and Table 4 (down-regulated). Notably, several genes related to chondrogenesis, bone repair, and differentiation were profoundly increased in response to combination treatment when compared to individual treatments. As show in the visualization of just the 1 h time points (Figure 3B), the co-expressed genes in patterns 8, 9 and pattern 11 vary depending on whether topotecan was given before or after the vehicle or oxaliplatin. On other hand, the coexpressed genes in patterns 14-16, vary depending on whether oxaliplatin was given before or after the vehicle. The 244 probes in patterns 8, 9, and 11, representing 59 co-expressed genes, enrich for pathways related to mRNA splicing, splicesomes, metabolism, cell cycle and DNA replication. The 188 probes in patterns 14-16, representing 45 co-expressed genes, enrich for pathways related to chromosome organization, chromatin packaging and remodeling. Full lists of genes derived from each pattern with absolute fold change of 4 or greater from their respective controls and p < 0.05, are provided in Supplementary Table 1.

PATHWAY ANALYSIS OF GENES WITHIN PATTERNS EXTRACTED BY EPIG
Gene lists within patterns revealing clear differences between control and treatment groups for each time-point were used to obtain pathway enrichment profiles by gene set enrichment analysis (Subramanian et al., 2005). Supplementary Tables 2-4 are lists (1 h, 6 h and 24 h respectively) of pathways from gene set enrichment analyses carried out against reference pathways, provided by Ariadne, within the designated patterns. Pathways common to multiple patterns were also identified and shown in Figure 4. The pathways commonly enriched in all patterns evaluated at all treatment times were those pathways related to chromatin remodeling and cell cycle regulation. Pathways derived from the lists of genes obtained from bone marrow samples 1 and 6 h after dosing were particularly enriched with pathways related to DNA repair, histone biology, cell cycle regulation, hypoxia, glutathione metabolism, and respiratory and oxidative phosphorylation. These regulatory events provide evidence of target-mediated biology for the drug treatments that can be potentially used as a basis for additional toxicodynamic modeling as early as 1 h after administration of a single dose.

DIFFERENTIALLY EXPRESSED GENES
Using statistical comparisons of treatments to time-matched controls (absolute fold change >2 and FDR < 0.01), 3304 gene probes in total were detected as differentially expressed (Supplementary Table 5). Principal component (PC) analysis of the gene expression data from the differentially expressed gene probes projected the samples in 3-dimensional space (Supplementary Figure 5) and revealed that PC #1 separates the samples by time and the top 3 PCs grouped the biological replicates by treatment very well.

ENRICHMENT OF BIOLOGICAL PROCESSES
Genes detected as differentially expressed at each time point (Supplementary Table 5) were used to enrich Gene Ontology (GO) biological processes (BPs). There were 641 GO BPs (Supplementary Table 6 Oxali is oxaliplatin (15 mg/kg) and Topo = topotecan (4.7 mg/kg). Veh1 is the vehicle used for Oxali, and Veh2 is the vehicle used for Topo. The data is the log base 2 ratio (treated sample to the average of the time-matched control) and the scale on the bottom displays the color range for the log base 2 ratio values. Red denotes upregulation, blue downregulation, and gray relatively no change.
directed acyclic graph derived from the terms in the cluster. DNA damage-signal transduction by p53 was highly enriched by the oxaliplatin exposure at the 6 h time point when it was given first but not second. Oxaliplatin given second elicited an ATP catabolic process at 6 h and positive regulation of epithelial to mesenchymal transition at 24 h. Topotecan when given second at 6 h impacted Ras GTPase activity in a positive regulation manner. The GTP catabolic process was enriched at the 6 h time point regardless of the order of administration of topotecan or oxaliplatin. Very few BPs were enriched at the 1 h time points. However, the biological processes highly connected to ventricular cardiac muscle cell development were enriched in a time-dependent manner, maximizing at the 24 h time point.

REGULATION OF p53 SIGNALING, APOPTOSIS AND CELL CYCLE RELATED GENES
Pathways related to DNA damage and p53-mediated cell cycle arrest were identified as highly connected and were uniquely Frontiers in Genetics | Toxicogenomics February 2015 | Volume 6 | Article 14 | 8 enriched in some samples. We then visualized the changes in gene expression within the rat p53 signaling Wikipathway at the 24 h time point in the study (Figure 6). As indicated by the relative expression to time-matched controls, MDM2, GADD45, SCOTIN, CASP8, and IGFB3 were up-regulated whereas the cyclins, p53, SIAH, GTSE1, and PARP1 were down-regulated. Regulation of several genes that also play roles in the apoptosis pathways (Figure 4) and cell cycle regulation (Figure 4) was also observed at the 24 h time point.

p53 SIGNALING PATHWAY GENE REGULATORY NETWORK RECONSTRUCTION
The enrichment of GO biological processes from the samples where oxaliplatin is given first followed by vehicle for topotecan revealed DNA damage, regulation of p53 signaling transduction (Figure 5). In addition, the overlay of gene expression data on the p53 signaling pathway revealed regulation of key components of the cascade at 24 h (Figure 6). We therefore sort out to reconstruct the gene regulatory network based on the 19 genes mapped to the rat p53 signaling Wikipathway. This would allow us to compare the gene interactions from the time point data in the samples where topotecan is given first followed by oxaliplatin vs. when oxaliplatin is given first followed by topotecan. Using the Gene Regulatory Network Inference (GRNInfer) software with the default setting to control the sparseness and the complexity of the network reconstruction, gene networks based on the average of the four replicate time point studies for each order of administration were revealed (Figure 7). When topotecan is given first followed by oxaliplatin, MDM2 proto-oncogene, E3 ubiquitin protein ligase (MDM2) and GADD45g are central hubs interacting with p53, cyclin-dependent kinases, several cysteine-aspartic acid proteases (CASPs), BID, two CASPs (CASP8 and CASP3) and other components ( Figure 7A). On the other hand, when oxaliplatin is given first followed by topotecan, CASP8 and the G-2 and S-phase expressed 1 gene (GTSE1) are the central hubs of the network interacting with p53, the cyclin-dependent kinases, CASP9, Kras, FAS, BID and other components ( Figure 7B). Essentially the activation and inactivation shown for components in the networks are caused by different central regulators depending on the order of administration.

DISCUSSION
This study was conducted to determine whether a single dose rodent study and toxicogenomics profiling informs the potential risk for enhanced bone marrow toxicity for the combined administration of clinically effective chemotherapeutic agents. Genomic profiling is a mature technology that has been strategically and efficiently used in preclinical drug safety assessment to predict safety issues that may be revealed in more lengthy, longer-term studies (Ryan et al., 2008;Huang et al., 2010). In addition, global gene expression data have shown promise in generating hypotheses about early onset "toxicity triggers" (Hamadeh et al., 2010).
Topotecan is a particularly difficult drug to use in combination because of the need for dose attenuation when given in combination regimens. Having a method to identify the potential for enhanced bone marrow toxicity when topotecan is given in combination with a second agent would be particularly useful.  (1), i.e., p = 1), Clustering of the pathways were performed using Pearson correlation dissimilarity (r) and average linkage grouping. Clusters (those with a 1-r ≥ 0.9) of GO BP terms (n ≥ 30) were labeled according to the node having the maximum number of paths to it within the GO BP subtree directed acyclic graph derived from the terms in the cluster. The color in the legend denotes the significance of enrichment. The more red the heat map color, the more significant the enrichment.
During the study, no drug-related adverse clinical signs were observed for any animals treated with topotecan or oxaliplatin alone or in combination. Histopathological evaluation of bone marrow established the presence of drug-related lesions in bone marrow (hypocellularity) in animals administered topotecan or oxaliplatin. Although our preclinical results highlight the need to consider the sequence of administration in a clinical protocol, our histopathological findings were not consistent with the sequence dependency of the severity of bone marrow toxicity that was reported in clinical studies of cisplatin and topotecan (De Jonge et al., 2000). In clinical trials, both neutropenia and thrombocytopenia were more severe when cisplatin was administered prior to topotecan and there were significantly lower absolute neutrophil count nadirs and percentage decrements in neutrophil and platelet counts with this sequence. Our results are consistent with the observation that oxaliplatin exhibits a favorable toxicity profile with a substantially lower incidence and severity of nephrotoxicity, ototoxicity, and myelosuppression.
Topotecan and other anticancer agents may cause DNA damage via mechanisms other than direct binding to DNA replication machinery. For example, camptothecin and analogs have been reported to induce apoptosis through a mechanism involving reactive oxygen species and oxidative stress pathways (Li et al., 2009). Similarly, pretreatment of mice with the anti-oxidant quercetin was reported to reduce topotecan-induced genotoxicity and cytotoxicity (Bakheet, 2011). Quercetin also reduced oxidative stress markers in topotecan treated bone marrow cells (Bakheet, 2011). Thus, it is feasible that the early effects of topotecan we observed in bone marrow are related to oxidative stress. Although it is clear that most of the in vivo biology for topotecan would be expected to be mediated by TOP1, the observed enrichment of genes related to ROS neutrophil-mediated cell damage and hypoxia-induced mitochondrial damage pathways, support additional mechanistic possibilities.
Traditional gene expression analysis has used two samples statistical tests for each comparison of interest or analysis of variance (ANOVA) model when the data are from studies with a factorial design. Upon selecting differentially expressed genes (DEGs), the next step is usually pooling the lists of DEGs to cluster analyze the data for visualization of patterns of gene expression across the samples. A limitation to this strategy is that the process of selecting the genes before detecting the patterns can omit genes that have salient expression profiles correlated across the samples. This caveat is of major concern when analyzing expression data comprised of genes whose expression are altered at low drug exposures or perturbed only for a short duration. The EPIG approach  was developed specifically for toxicogenomics and other series type studies where the algorithm uses an ANOVAlike statistical evaluation of the patterns of gene expression but also takes into account correlation of genes within an expression pattern. In the EPIG two-step approach, all significant patterns in the data are extracted first, followed by categorization of gene expression profiles. Application of the EPIG method to our microarray data and subsequent pathway and gene regulatory network analyses allowed identification of genes, pathways and regulatory interactions that appear to represent promising biomarker signatures for bonemarrow toxicity. Pathways related to p53 signaling, DNA repair, histone biology, cell cycle regulation, hypoxia, glutathione metabolism, and respiratory and oxidative phosphorylation reflect the biological effects of either treatment on the bone marrow. The considerable enrichment of pathways involved with DNA repair and chromatin remodeling is remarkably wellaligned with the biological mechanisms and downstream effects of topotecan. The interplay between chromatin remodeling and DNA repair factors is infrequently discussed in relation to DNA damage response mechanisms of the bone marrow. Our analysis highlights a meaningful relationship between chromatin remodeling complexes and mechanisms of bone marrow toxicity and repair that warrants further investigation.
An interesting biological response pathway of genes related to tissue injury was derived from our analysis of co-expressed genes (Supplementary Table 1). Regulation of chondrogenesis, bone repair, and differentiation was identified for some of the genes that were profoundly increased in response to combination treatment when compared to individual treatments. These genes were matrix metallopeptidase 12 (MMP12); transgelin (TAGLN, SM22α); cyclin D1 (CCND1); serpin peptidase inhibitor, clade H (heat shock protein 47), member 1 (collagen binding protein1) (SERPINH1); tenascin C (TNC); ADAM metallopeptidase with thrombospondin type 1 motif, 9 (ADAMTS9); and bone morphogenetic protein 2 (BMP2). TNC is an extracellular matrix glycoprotein that is specifically and transiently expressed upon tissue injury. SERPINH1 functions as a molecular chaperone during collagen synthesis and maturation (Nagata, 1998;Lamande and Bateman, 1999;Razzaque and Taguchi, 1999;Hendershot and Bulleid, 2000). Upon tissue damage, TNC regulates a wide variety of pathways that mediate both inflammatory and fibrotic processes, enabling effective tissue repair (Truong et al., 1996;Chiquet-Ehrismann and Chiquet, 2003;Midwood and Orend, 2009). TAGLN is a shape change sensitive 22 kDa actin-binding protein of the calponin family that may regulate conversion of FIGURE 7 | Reconstruction of the p53 signaling gene regulatory network. The 19 genes mapped in the rat p53 signaling Wikipathway is used to reconstruct the network based interactions derived from the gene expression data. For (A) topotecan followed by oxaliplatin and (B) oxaliplatin followed by topotecan, the average of the biological replicates' gene expression data at each time point were used to reconstruct the gene regulatory network depending of the order of administration of the two drugs. A red arrow indicates activation, a blue arrow indicates inactivation. The labeling of the nodes is based on the UniGene symbol.
adult bone marrow-derived mesenchymal stem cells into smooth muscle cells (SMCs) and is an early marker of smooth muscle differentiation (Lawson et al., 1997). Expression of ADAMTS9 was shown to be up-regulated during chondrogenic differentiation of human mesenchymal stem cells (Boeuf et al., 2012).
Gene profiling of bone marrow environment cells revealed distinct expression profiles for genes encoding for ADAMs and their inhibitors (Bret et al., 2011). In the current study, all cells expressed ADAMTSs genes at a low level, with the exception of bone marrow stromal cells. BMP2 plays an essential role in chondrocyte proliferation and maturation during endochondral bone development (Shu et al., 2011). Similarly, MMPs are required for both endochondral and intramembranous ossification during bone repair and it is likely that this gene response is a sensitive indicator of initial degradation of extracellular matrix. These gene responses are similar to microarray studies that identified gene responses during stages of bone marrow-ablation-induced bone regeneration (Wise et al., 2010). Taken together, the observed gene responses may represent a unique biomarker panel for bone marrow that will flag early tissue damage onset followed by chondrogenesis and intramembranous regeneration processes. The 11-and 7-fold increases in MMP12 and TAGLN observed 24 h post-treatment with topotecan compared with a 54-and 17-fold increase when topotecan treatment is given with oxaliplatin might represent an additive response for the combination. The effect of treatment combinations on these responses likely reflects a change in the extent of damage and response of the bone marrow to injury.
Among the genes showing the most profound decreases in response to all treatments are several genes for proteins regulating shape and hemolysis of erythrocytes. We noted a 26-fold decrease in spectrin, alpha, erythrocytic 1 (elliptocytosis 2) (SPTA1). Mutations in spectrin genes that render red cells deficient in spectrin are associated with abnormal cytoskeletal architecture making erythrocytes susceptible to hemolysis (Kakhniashvili, 2001;Broderick and Winder, 2005). A similar contribution to actin dynamics in platelets has been elucidated for pleckstrin-2 (PLEK2) (Lian et al., 2009), and we report a 14-fold decrease in PLEK2 mRNA in response to topotecan alone. In addition, a 9fold decrease was noted in Ankyrin 1 (ANK 1), an erythrocyte membrane protein that is defective in many patients with hereditary spherocytosis, a common hemolytic anemia. Taken together our results demonstrate a clear connectivity between the most profoundly affected genes and the clinical adverse effect profile for topotecan. For example, a recent clinical study designed to determine the dose of weekly oral topotecan allowing safe administration in patients with recurrent gynecologic malignancies, 13 (11.1%) doses of drug were held because of anemia in 8 patients, neutropenia in 7, or thrombocytopenia in 2 (Von Gruenigen et al., 2012).
This study reports gene response data in bone marrow posttreatment with topotecan or oxaliplatin alone and in combination following a single administration. There are no comparable gene expression studies on the effect of the combination of the two agents used here, but one study, based on 8 mg/kg oxaliplatin in rat bone marrow at 24 h, was found in a public data source. That public (yet unpublished) study is deposited in the Drug Matrix www.frontiersin.org February 2015 | Volume 6 | Article 14 | 17 8.0 database (Ganter et al., 2005) hosted by the NIEHS/NTP at https://ntp.niehs.nih.gov/drugmatrix/index.html. One published study compares the effect of oxaliplatin in ovarian cancer spheroids (L'Esperance et al., 2008). This underscores the uniqueness and novel aspect of our combination study. Supplementary  Figures 6, 7 show the paucity of overlap of the oxaliplatin DEGs from our current study and those from the aforementioned public studies (Drug Matrix bone marrow and ovarian cancer spheroids respectively). Caution is needed in interpreting these comparisons due to the differences in the doses administered, the target cell type/tissue, the array platforms, significance test and the annotations based on gene symbol to compare results. We provide a comprehensive list of genes from our microarray analysis of bone marrow taken at various time points after treatment of rats to enable additional analysis and hypothesis generation (http:// tools.niehs.nih.gov/cebs3/ui/). Although we will continue to analyze these data for biological response pathways and mechanistic pathways for toxicity, we publish these data to encourage other investigators to employ alternative analysis methods and propose additional relationships relevant to the observed bone marrow responses and treatments. The use of a single dose study using genomic endpoints as a measure of enhanced bone marrow toxicity is exciting for several reasons. This approach could be incorporated into a mechanismbased risk evaluation and a single dose administration can be readily incorporated as an in vivo or in vitro screening paradigm. For example, attempts to validate in vitro bone marrow systems could use these data to explore and identify molecular anchors that translate from in vitro to in vivo. Second, this proof of concept revealed that the molecular responses in bone marrow toxicity is much earlier than previously documented preclinical histopathological observations of topotecan toxicity. Examining the acute response of a target organ to chemotherapy and obtaining a molecular footprint associated with this response in combination with a targeted agent can guide us along new avenues of investigation about mechanisms of target-mediated toxicity. At this point, these data represent qualitative information that is relevant to understanding mechanisms of toxicity. The initial phenotypic anchoring of these data with additional endpoints may ultimately provide a more comprehensive understanding of pathways underlying bone marrow toxicities. Indeed, other reports have highlighted the potential use of toxicogenomics data to enable integrative risk assessment and biomarker identification (Ellinger-Ziegelbauer et al., 2011;Matheis et al., 2011). We plan to apply this approach to other target organs and additional anticancer drug combinations to examine the broader implications of this strategy.