Genomic analysis of microRNA time-course expression in liver of mice treated with genotoxic carcinogen N-ethyl-N-nitrosourea

Background Dysregulated expression of microRNAs (miRNAs) has been previously observed in human cancer tissues and shown promise in defining tumor status. However, there is little information as to if or when expression changes of miRNAs occur in normal tissues after carcinogen exposure. Results To explore the possible time-course changes of miRNA expression induced by a carcinogen, we treated mice with one dose of 120 mg/kg N-ethyl-N-nitrosourea (ENU), a model genotoxic carcinogen, and vehicle control. The miRNA expression profiles were assessed in the mouse livers in a time-course design. miRNAs were isolated from the livers at days 1, 3, 7, 15, 30 and 120 after the treatment and their expression was determined using a miRNA PCR Array. Principal component analysis of the miRNA expression profiles showed that miRNA expression at post-treatment days (PTDs) 7 and 15 were different from those at the other time points and the control. The number of differentially expressed miRNAs (DEMs) changed over time (3, 5, 14, 32, 5 and 5 at PTDs 1, 3, 7, 15, 30 and 120, respectively). The magnitude of the expression change varied with time with the highest changes at PTDs 7 or 15 for most of the DEMs. In silico functional analysis of the DEMs at PTDs 7 and 15 indicated that the major functions of these ENU-induced DEMs were associated with DNA damage, DNA repair, apoptosis and other processes related to carcinogenesis. Conclusion Our results showed that many miRNAs changed their expression to respond the exposure of the genotoxic carcinogen ENU and the number and magnitude of the changes were highest at PTDs 7 to 15. Thus, one to two weeks after the exposure is the best time for miRNA expression sampling.


Background
MicroRNA (miRNA) is a class of small nucleic acids that negatively regulate gene expression [1]. They are singlestranded RNA molecules that are not translated into proteins. Each miRNA molecule is partially complementary to one or more mRNA transcripts, and functions to down-regulate gene expression by inhibiting protein translation or destabilizing target transcripts [2]. The expression of miRNAs is regulated developmentally and spatially, and is involved in differentiation and proliferation of cells [3]. Therefore, miRNA molecules can modulate a wide array of growth and differentiation processes in cancer [4]. A number of studies have demonstrated that miRNA expression is commonly dysregulated in human cancer; and that miRNAs are extensively involved in carcinogenesis and act as either dominant or recessive cancer genes [5]. miRNA expression profiling has shown promise in defining tumor malignancy status and is surprisingly informative when used to identify tumor types, differentiation states and developmental lineages [6].
The available information on miRNA function suggests that miRNA expression profiles might also have predictive value for assessing chemical carcinogenicity. It has been reported that expression of some miRNAs is associated with tumor initiation [7]. Thus, specific miRNAs could represent attractive molecules as informative biomarkers of exposure to carcinogens. Studies on the relationship between miRNAs and carcinogen exposure have been previously reported [8] and miRNA expression has been shown to be dysregulated by many carcinogenic agents like 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanone (NNK) [9], 7,12-dimethylbenz[a]anthrance [10], 2-acetylaminofluorene (2-AAF) [11] and radiation [12]. The results of these studies indicate that miRNAs are involved in early stages of carcinogenesis and suggest that miRNAs could be a useful tool for detecting carcinogen exposure. Expression data from these studies, however, were generated primarily from chronic or subchronic carcinogen exposures in which the animals were treated with multiple doses of carcinogens over a long period such as 20 or 24 weeks [9,11]. This kind of studies are able to reveal cumulative effects of chemical treatments on miRNA expression. However, they are unable to provide time course miRNA expression data caused by a single dose of chemical treatment.
Time-course data is an important component of toxicological studies. It reveals a toxicological response as a transient, continuous, or delayed response. Since the regulation of miRNA expression is a dynamic process, a temporal design provides information regarding an appropriate sampling time for miRNA expression change after carcinogen treatment. In this study, N-ethyl-Nnitrosourea (ENU), a model genotoxic carcinogen, is used for the carcinogen treatment. We treated mice with one dose of ENU and measured the expression level of miRNAs in the liver of the treated and control mice at several posttreatment times. This treatment and sampling design allowed us to analyze the time-course changes of miRNA expression in tissues exposed to the carcinogen. We have previously published the mutant frequency data from mice treated in this study [13] ENU is mutagenic in a wide variety of mutagenicity test systems and carcinogenic in various organs of mammals [14]. It induces hepatocellular carcinomas in mouse liver that receives a single dose of ENU [15]. Thus, ENU is a suitable genotoxic carcinogen for studying the time course expression of miRNAs. Usually, biological effects of a carcinogen are affected by various factors, including chemical absorption, distribution, metabolism, and elimination. ENU, however, directly alkylates nucleotides without metabolic activation [14]. Thus, choosing ENU as a model carcinogen avoids the effects of these unrelated processes to a large extent due to ENU's direct activity. Most importantly, abundant information on the time-course effects of ENU toxicity and carcinogenicity has been accumulated including data regarding DNA adducts [16], gene mutations [13], gene expression [17] and tumor formation [18]. This information is valuable to understand alteration of miRNA expression induced by ENU exposure.

Principal component analysis
To globally view the temporal changes of miRNA expression induced by a single ENU treatment at different posttreatment time points, principal component analysis was performed to classify all of the 35 animal samples (5 samples each group for post-treatment days [PTDs] 1, 3, 7 and 15 after ENU treatment; 4 samples each group for PTDs 30 and 120; and 7 samples for the control group). The normalized Ct values of 376 miR-NAs (Additional File 1) were used for this analysis and the results for the first 3 principal components are shown in Figure 1. The samples can be roughly divided into 2 groups. Samples at PTDs 7 and 15 are grouped together while samples from the other groups including the control group and treatment groups at PTDs 1, 3, 30 and 120 are clustered together. The results from the principal component analysis suggest that the miRNA expression was globally altered by ENU treatment at PTDs 7 and 15, but not at PTDs 1, 3, 30 and 120 although two samples at PTD 3 also fall into the PTDs 7 and 15 group, which might reflect that some animals at PTD 3 had miRNA expression patterns similar to those from samples at PTDs 7 and 15.

Time course analysis of miRNA expression
There were 43 miRNAs whose expression was dysregulated by ENU treatment at least at one posttreatment time. The names of those differentially expressed miR-NAs (DEMs), their fold changes against the controls and the posttreatment time points when their expression was altered, as well as their possible functions reported in the literature are listed in Table 1 (see Additional File 2 for the fold changes and p values of all the miRNAs in the PCR array). The amount of the DEMs was changed with time after ENU treatment ( Figure 2). Only 3 and 5 DEMs were found at PTDs 1 and 3, respectively. This amount increased to 14 and 32 at PTDs 7 and 15, and then decreased to 5 each at PTDs 30 and 120. Most of the DEMs (4 out of 6) at PTDs 1 and 3 were also observed at PTDs 7 and 15. Two miRNAs, mmu-miR-34a and mmu-miR-34b-5p, were up-regulated at all posttreatment time points except day 120. Most of the miRNAs that were differentially expressed at both PTDs 7 and 15 had the same change directions and similar change magnitudes at the two time points, indicating miRNA alteration during this period was relatively stable (Table 1). However, DEMs at PTD 120 were different from those found in previous posttreatment time points except for mmu-miR-453 that was also dysregulated at PTDs 7 and day 15 but with an opposite change direction ( Table 1).
Confirmation of the temporal expression changes of three miR-34 family miRNAs and one miR-762 family miRNA by individual TaqMan assays The miR-34 family is worth special attention because two of its members, mmu-miR-34a and mmu-miR-34b-5p, were significantly up-regulated one day after ENU exposure and maintained increased expression at the 5 subsequent time points up to PTD 30, while another family member, mmu-miR-34c, displayed significant over-expression at multiple time points from PTD 3 to 30. To serve as a confirmation of the data from RT 2mouse miRNA PCR Array, a different real-time PCR platform, TaqMan miRNA assay, was used to measure the time-course changes of miR-34 family members induced by the ENU treatment using the same RNA samples. For TaqMan miRNA assay, all of the samples were reversely transcribed simultaneously and the PCRs for the different samples were run in the same plates. This procedure can greatly reduce the inter-plate variability that could be produced when PCR arrays are used. A comparison of miR-34 family miRNA expression measured by the two different platforms is shown in Figure  3. The results from the two real-time PCR assay platforms are very consistent and show similar temporal kinetics of miRNA expression for miR-34 family miR-NAs, rising from day 1 or 3, reaching peaks at day 15, and decreasing until the end of observation, day 120. Another miRNA, mmu-miR-762 that is not similar with miR-34 family miRNAs in sequence, were also examined to confirm the array data. The consistency between two platforms was also observed ( Figure 3).

Hierarchical clustering analysis of the DEMs at different posttreatment times
The 43 DEMs were grouped using hierarchical clustering methods based on their logarithmically transformed fold changes. The result is shown in Figure 4. The DEMs were clustered into four groups. Groups I and II included the DEMs whose expression was mainly upregulated. While the DEMs in Group II were up-regulated across all time points, those in Group I showed slightly down-regulated in one or two time points. Groups III and IV composed of DEMs whose expressions were down-regulated. Day 7 and 15 clearly show the highest fold changes compared to other time points, indicating miRNA expression reached maximal levels or changes between the two time points (Figure 4 and 5).

In silico Functional analysis of the DEMs
A literature search of the PubMed database was conducted using the gene name of each DEM without considering species difference since miRNA functions are relatively conservative across different species [19]. Twenty seven out of 43 of the DEMs have been functionally investigated by biological experiments ( Table 1). Many of these miRNAs were dysregulated in tumors or  let-7b + Ectopic expression of let-7b reduced HMGA2 expression and cell proliferation in a lung cancer cell line [33].
miR-106b + miR-106b was up-regulated in several human tumors compared with adjacent normal tissues and formed a negative-feedback loop with cell cycle regulator E2F1 [43].
miR-130a + Expression of miR-130a was significantly up-regulated in primary glioblastomas compared with normal peripheral brain tissue [44].
miR-487b + ND miR-590-3p + + + miR-590-3p and other miRNAs were suggested to mediate control of autoimmune gene expression [61]. miR-590-5p + miR-590 was involved in regulating the expression of transforming growth factor TGF-beta1 and its receptor TGF-betaRII [62]. miR-672 + + + ND miR-677 + + ND miR-700 + ND tissues exposed to carcinogens. They played roles in cell proliferation, cell cycle regulation, cellular transformation, immune response, invasion, apoptosis, tissue development and growth. Some of the DEMs are putative tumor suppressors or oncogenes such as miR-34a (see literature cited in Table 1). Although the functions of the DEMs are highly related to tumorigenesis according to the reports, the studies on the miRNA functions were conducted in different species and various biological systems. To explore the functional relationship between miRNA expression changes and ENU treatment, we used the top computationally-predicted target genes of these DEMs according to the sequence complementarity between miRNA and their target mRNAs for our functional analysis. Such approaches have been widely used for functional analysis of miRNA target genes since miRNAs exert their gene regulatory activity primarily by base pairing to the 3' UTR of their target mRNA, leading to mRNA degradation or translational inhibition [20]. Thirty six DEMs at PTDs 7 and 15 were used for this analysis because the principal component analysis showed that the global miRNA expression was altered by ENU treatment at these two time points. The selected target genes of the miR-142-3p + miR-142-3p was over-expressed in childhood B-cell precursor acute lymphoblastic leukemia [64].
Note: miRNAs with p value < 0.01 and fold change >2.0 were considered significantly changed. The symbols +, ++, and +++; or -, --, and ---indicate that the miRNA was up-or down-regulated at a magnitude of 2-4, 4-6, and >6 fold, respectively. No symbol is given if the miRNA was not significantly changed. ND means not determined.   DEMs were input into IPA, and the related functions and pathways were determined and are provided in Figure 6 and Additional File 3. The most relevant biological functions putatively controlled by the ENU-induced DEMs are related to mutagenesis and tumorigenesis, involved in tumor/cancer formation, cell cycle and DNA repair, as well as other ENU-toxicity related functions.

Discussion
To explore the temporal response of miRNA to treatment of genotoxic carcinogens, mice were administrated with one dose of 120 mg/kg ENU, which can significantly induce tumors and mutations in mouse liver [13]. The ENU treatment resulted in temporal changes in miRNA expression in mouse liver and the altered  Temporal changes in miRNA expression after ENU treatment miRNA aberrant expression has been found to be a common feature of tumor malignancy [21]. Several studies on carcinogen-altered miRNA expression have been conducted and showed that chronic or subchronic carcinogen exposure can dysregulate miRNA expression. Long-term exposure of female Fisher 344 (F344) rats to a tamoxifen-containing diet led to alterations in the miRNA expression profile in liver tissue prior to tumor formation [22]. Sprague-Dawley rats treated with 2-AAF for 12 or 24 weeks exhibited disrupted regulation of the miR-34a-p53 feed-back loop and substantial deregulation of expression of miR-18, miR-21, miR-182, and miR-200 family miRNAs [11]. Male F344 rats continuously fed with NNK for up to 20 weeks resulted in alteration of miRNA expression in the lungs of rats [9]. Although these studies demonstrated that carcinogen exposures were able to alter expression of certain miR-NAs, such experiments could not provide information on temporal response of miRNAs to carcinogen exposure. The alteration of miRNA expression by these subchronic to chronic treatments do not provide the response period of miRNAs to treatment. The changes could result from the cumulative effects of persistent long-term chemical exposure or from a short term response. At present, the duration of miRNA response to exposure to a genotoxic carcinogen in animals or humans is unknown. In this study, mice received only a single dose of a potent genotoxic carcinogen ENU and were sacrificed at different time points so that temporal response of miRNA expression to the treatment can be determined. The results clearly demonstrate that exposure of mice to a single dose of ENU can cause distinctive alterations in miRNA expression at the different sampling times after the ENU treatment (Table 1 and Figures 1, 2, 3 and 4). Sampling time after the ENU exposure appears to be a very important parameter in miRNA expression in both the amount of DEMs and their fold-change levels. Principal component analysis indicates that mouse livers sampled at 7 and 15 days after ENU treatment display distinctive miRNA expression patterns compared to controls and other sampling time points (Figure 1). Three DEMs were found one day after the ENU treatment. The number of DEMs increased with the sampling time, reaching the peak at PTD 15 with 32 DEMs. The amount of DEMs declined after 15 days and only 5 DEMs were identified at PTDs 30 and 120 ( Figure 2). The magnitude of the DEMs' expression alteration was also changed with time. Most of the DEMs had their highest fold-changes at PTDs 7 or 15 (Figures 4 and 5). Hierarchical clustering analysis of the DEMs divided them into four groups according to the change direction and intensity of these DEMs' expression alteration. The DEMs in Groups II and IV responded to the ENU exposure quickly, indicating that some miRNAs can be changed by the treatment within a few days. However, regarding the alteration of miRNA expression induced by the ENU treatment, both the amount and the change magnitude of the DEMs peaked at PTD 15.
The pattern of temporal changes in miRNA expression after ENU treatment is different from previous reports on gene expression changed by ENU treatment [17]. In the previous studies, gene expression changes in both amount and intensity peaked at 4 hours and declined at 20 hours, 14 and 28 days after administration of ENU. The mechanism(s) for the difference between the temporal changes induced by ENU in miRNA expression and gene expression is unknown. It is possible that many genes, such as DNA repair genes, can directly respond to the DNA damage caused by ENU exposure whereas miRNAs as posttranscriptional regulators for gene expression may respond indirectly to ENU insults via the alteration of gene expression [23]. Also, turnover rates for mature mRNAs and miRNAs could play a role in the difference. Mature miRNAs generally have lower turnover rates and exist longer than mRNAs [24].
The pattern of the temporal expression change could result from biological responses of liver miRNAs to the genotoxic and cytotoxic effects of ENU. miRNA expression may vary between different biological mechanisms like DNA adduct formation by ENU, cell death induced by the DNA damage, and cell proliferation due to the cell death. The DNA damage induced by ENU is very fast, occurring within hours after treatment [25]. However, the DNA repair and cell proliferation processes could take a few weeks [26]. Some miRNAs like miR-34a could respond to the DNA damage quickly and exhibited dysregulation within one day while most miR-NAs that function mainly as regulators for cell differentiation and proliferation altered their expression only after several days.
In silico pathway analysis of the functions of the ENUinduced DEMs ENU, as a potent monofunctional ethylating agent, reacts directly with the nucleophilic nitrogen and oxygen atoms in DNA and with the oxygen atoms in the backbone phosphates, forming various ethylated products including N 7 -ethylgunine, N 3 -ethyladenine, O 6ethylguanine, O 2 -and O 4 -ethylthymine, and O 2 -ethylcytosine [16,25]. Once formed in the DNA, these products can result in cell death, like apoptosis [27] or become the substrate of cellular repair processes of various kinds, such as specific dealkylation by the O 6 -alkylguanine-DNA-alkyltransferase (AGT) [28], removal of N 7alkylguanine or N 3 -alkylguanine by specific glycosylases [29], or removal O 4 -alkylthymine by the more general action of the nucleotide excision repair system [30]. Also, ENU is known to induce cell proliferation [31]. If miRNAs are involved in the regulation of the genes that are related to genotoxic functions like cell proliferation, cell cycle arrest, apoptosis, and DNA repair, expression of miRNAs related to these functions should be changed by the ENU treatment. Indeed, our literature search results indicate that most of the DEMs induced by ENU treatment function as regulators for cell cycle arrest [32], cell proliferation [33], apoptosis [34], DNA repair [35] and other biological processes related to ENU cytotoxicity, genotoxicity and carcinogenicity (Table 1).
Because each miRNA can regulate many target genes and several miRNAs may affect a single gene, it is important to analyze the functions of all DEMs together. Computational approaches have been a major focus in determining the general principles that are thought to govern miRNA target recognition and mode of action. In this study, the target genes computationally predicted by the miRanda algorithm were used for functional analysis. The miRanda algorithm was developed by the Sanger Institute and is widely used in miRNA studies [36]. The top 5% target genes of all of the DEMs were selected and used for IPA functional analysis. The top functions affected by the ENU exposure are tumor morphology, cell cycle, DNA replication, recombination and repair, and cancer. These biological functions identified by the analysis show that the DEMs are related to ENUinduced carcinogenesis ( Figure 6). For example, many genes involved in DNA repair or response to DNA damage can be dysregulated by ENU treatment [37] and the miRNAs that target these genes could change their expression to regulate these functional processes.

MiR-34 family might have the potential to be explored as a biomarker for genotoxin exposure
Our results indicate that some miRNAs responded to ENU treatment with a wide temporal range. These miR-NAs might have the potential to be used as biomarkers for predicting the genotoxic carcinogenicity of chemicals. Among these miRNAs, the miR-34 family is worth special attention. All of the 3 miRNAs were significantly changed at four different time points (Figure 3). Their expressions were enhanced by 3.21-fold (miR-34a), 3.11fold (miR-34b) and 2.37-fold (miR-34c) on PTD 1 and the fold changes continued to increase and peaked at PTDs 7 or 15. The miR34 family genes are the direct transcription targets of tumor suppressor p53 [32,38]. miR-34b and miR-34c are encoded by the same primary transcript from chromosome 11 in human or chromosome 9 in mouse while miR-34a is located in a different chromosome [35]. The promoter region of miR-34a and miR-34b/c each contain a palindromic sequence that matches the canonical p53 binding site and can be bound by p53 as shown by chromatin immunoprecipitation [32]. Interestingly, our results found that miR-34b and miR-34c changed in correlated manner at all the sampling time points (Figure 3). miRNAs in miR-34 family play important roles in various p53-initiated biological processes. Up-regulation of miR-34a and miR-34b/c caused a cell-cycle arrest in the G1 phase [32]. miR-34b/c inhibits cell proliferation and colony formation in soft agar [39]. Introduction of miR-34a and miR-34b/c into primary human diploid fibroblasts induces cellular senescence [35]. Re-expression of miR-34a in tumor cells induced apoptosis [38]. These biological processes controlled by miRNAs in the miR-34 family are related to ENU cytotoxicity, genotoxicity, and carcinogenicity. Our results indicate that the miR-34 family of miRNAs seems to have the potential to be valuable biomarkers for toxicological application.

Conclusions
Our study indicates that one dose treatment of ENU, a chemical inducing tumors and mutations, resulted in deregulation of a large number of miRNAs. In silico functional analysis suggested that these miRNAs were related to ENU mutagenesis and carcinogenesis in the mouse liver. The deregulation of ENU-induced miRNA expression changed with time and peaked at day 15 after the treatment. The findings suggest that one to two weeks after ENU exposure is the best time for miRNA expression sampling. Moreover, miRNAs in the miR-34 family worth further study to explore their potential as biomarkers for exposure of genotoxic carcinogens.

Animal treatment
The animal treatment protocol and mutant frequency analysis was described previously [13]. Briefly, sixmonth-old female Big Blue mice were injected intraperitoneally with a single dose of 120 mg/kg body weight ENU (CAS# 759-73-9, Sigma, St. Louis, MO) or the vehicle dimethylsulfoxide (DMSO, Sigma) in 1 ml/kg body weight (0.1%). For ENU treatment, groups of 4 or 5 animals were sacrificed on PTDs 1, 3, 7, 15, 30, and 120. For the vehicle treatment, 4 and 3 animals were sacrificed on PTDs 1 and 30, respectively. DMSO is not carcinogenic and mutagenic [40]. It has showed no effects on gene expression in the dose that we used [41]. Also, our previous study demonstrated that DMSO did not change mutant frequency in mouse liver at different sampling times [13]. Therefore, the 7 control samples were grouped together as a common control group for the treatment groups sampling at the different time points to increase the statistical power. The tissues were isolated and frozen at -80°C. The liver samples were used for this study. All animal experiments were conducted by following the recommendations set forth by our Institutional Animal Care and Use Committee.

miRNA isolation
About 60 mg of liver tissue was cut from each frozen liver sample and suspended in RNAlater-ICE (Ambion Inc., Austin, TX). The tissue pieces were transferred to 600 μl RNA lysis/binding buffer and minced using Tissue Tearor™ (BioSpec Products Inc., Bartlesville, OK). miRNAs were isolated using mirVana™ miRNA isolation kit (Ambion) that specifically captures small RNAs with length of less than 200 nucleotides. The isolated RNAs were resolved in 100 μl nuclease-free water (Ambion). RNA concentrations were determined using NanoDrop 1000 Spectrophotometer (NanoDrop Technologies, Wilmington, Delaware). The quality of RNA samples was characterized on an Agilent BioAnalyzer (Agilent Technologies, Santa Clara, CA) using an RNA6000 Nano Chip (Agilent).

PCR Array analysis of miRNA expression
Two hundred nanograms of enriched small RNA were converted into cDNA using RT 2 miRNA First Strand Kit (SABiosciences Corporation, Frederick, MD). The cDNAs were mixed with 2 × RT 2 SYBR Green PCR Master Mix (SABiosciences) and dispersed into 384-well Mouse Genome miRNA PCR Array (MAM-3100E, SABiosciences) with 10 μl/well reaction volume. The PCR array contained a panel of primer sets for 376 mouse miRNAs, four small RNAs as the internal controls and four quality controls. The real-time qRT-PCR was performed on a 7900 real-time PCR system (Applied Biosystems Inc., Foster, CA) with following cycling parameters: 95°C for 10 mins, then 40 cycles of 95°C for 15 s, 60°C for 30 s and 72°C for 30 s. SYBR Green fluorescence was recorded from every well during the annealing step of each cycle. The threshold cycle (Ct) value of each sample was calculated with software SDS 2.3 (Applied Biosystems). To calculate Cts, we set the threshold line as 0.15 and kept it the same across all of the analyses. The baseline was automatically defined by the software.

Normalization and statistical analysis
Normalization and statistical analysis of miRNA expression were conducted using SABiosciences Online PCR Array Data Analysis Web Portal. MiRNA expressions were compared between the treatment group at each time point and the control group. The ΔΔCt method was utilized to calculate the fold change (FC). Four genes, snoRNA251, snoRNA202, snoRNA142, and U6 in the PCR arrays, were averaged as the endogenous control and the vehicle control group was used as external control to normalize each sample. The formula: FC = 2[ -(mean of ΔCt values of treated samples -mean of ΔCt values of control samples)] was used for up-regulated gene, while FC = -2^(mean of ΔCt values of treated samples -mean of ΔCt values of control samples) was used for the down-regulations. T-tests were used to calculate the p value to determine whether there is a significant difference for miRNA expression between the control and the treatment groups for each miRNA at each time point. miRNAs with p < 0.01 and the absolute value of FC >2.0 were considered as DEMs.

Principal component analysis and hierarchical clustering analysis
Principal component analysis of expression profiles of all miRNAs from each time point after the ENU treatment was conducted using the autoscaled method within ArrayTrack [42]. The normalized ΔCt values were used for this analysis and the analysis was performed without filtering any miRNAs.
To examine types of expression changes, hierarchical clustering analysis was performed using R software http://www.r-project.org/. The miRNAs whose expressions were significantly differentially expressed at least at one sampling time were used for the clustering. The FCs for control samples were set to zero. Euclidean and Wards methods were used for distance-calculation and linkage, respectively.

In silico Functional analysis of the DEMs
The miRanda database was used for identification of DEMs' target genes http://microrna.sanger.ac.uk/ sequences/. The top 5% of the most reliable predicted target genes of the DEMs at PTDs 7 and 15 were selected according to the p values given in the database. A total of 1376 genes were determined as the predictive target genes of these DEMs. The selected target genes were then input into Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems, Inc., Redwood City, CA). Ingenuity Core Analysis and Knowledge Base was used as the reference set. IPA interpreted the genes in the context of biological processes, pathways and molecular networks and defined the most relevant biological functions of the predicted genes of the DEMs.
TaqMan qPCR confirmation of the temporal expression changes of miR-34 family miRNAs TaqMan MicroRNA Assays were used to confirm the temporal expression changes of 3 miR-34 family members, mmu-miR-34a, mmu-miR-34b-5p, and mmu-miR-34c, as well as a miR-762 family member, mmu-miR-762. The experiment also severed as the verification of the RT 2 -mouse miRNA PCR array assay. The TaqMan miRNA assay kits were purchased from Applied Biosystems (Foster City, CA). The same small RNA samples used for the PCR arrays were also used for the TaqMan miRNA assays. The experiment was performed by following the manufacturer's protocol. In brief, each 10 μl reverse transcription (RT) reaction contained 44 ng of small RNA, 50 nM stem-loop RT primer, 1 × RT buffer, 0.25 mM each of dNTPs, 3.33 U/μl MultiScribe™ reverse transcriptase and 0.25 U/μl RNase inhibitor. The RT reactions were incubated in a GeneAmp PCR System 9700 (Applied Biosystems) for 30 min at 16°C, 30 min at 42°C, followed by 5 min at 85°C, and then held at 4°C. Each real-time PCR reaction (10 μl volume) containing 0.78 μl of RT product, 5 μl of 2× TaqMan Universal PCR Master Mix, and 0.5 μl TaqMan MicroRNA assay (the mixture of TaqMan probe, forward primer, and reverse primer). The PCR reaction was conducted in an Applied Biosystems 7500 Fast Real-Time PCR System at 95°C for 10 min, followed by 40 cycles of 95°C for 15 sec and 60°C for 1 min. The threshold cycle (Ct) is defined as the fractional cycle number at which the fluorescence exceeds the fixed threshold of 0.02. Four samples at each time point were used for the TaqMan confirmation.

Accession numbers
All PCR array Ct raw data are available through Gene Expression Omnibus (Series accession numbers: GSE20248).

Additional material
Additional file 1: The normalized Ct values (dCt) for the control samples and treatment samples at different days after ENU treatment.
Additional file 2: The p values and fold changes of microRNAs at different days after ENU treatment.
Additional file 3: The biological functions related to the top 5% targets of DEMs.