Email updates

Keep up to date with the latest news and content from BMC Genomics and BioMed Central.

Open Access Research article

MicroRNA-mediated gene regulation plays a minor role in the transcriptomic plasticity of cold-acclimated Zebrafish brain tissue

Ruolin Yang1, Zhonghua Dai1, Shue Chen12 and Liangbiao Chen1*

Author Affiliations

1 State Key Laboratory of Molecular Developmental Biology, Institute of Genetics and Developmental Biology, Chinese Academy of Sciences, Beijing 100190, China

2 Graduate School of the Chinese Academy of Sciences, Beijing 100080, China

For all author emails, please log on.

BMC Genomics 2011, 12:605  doi:10.1186/1471-2164-12-605


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2164/12/605


Received:25 May 2011
Accepted:14 December 2011
Published:14 December 2011

© 2011 Yang et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

MicroRNAs (miRNAs) play important roles in regulating the expression of protein-coding genes by directing the degradation and/or repression of the translation of gene transcripts. Growing evidence shows that miRNAs are indispensable player in organismal development with its regulatory role in the growth and differentiation of cell lineages. However, the roles of miRNA-mediated regulation in environmental adaptation of organisms are largely unknown. To examine this potential regulatory capability, we characterized microRNAomes from the brain of zebrafish raised under normal (28°C) and cold-acclimated (10°C, 10 days) conditions using Solexa sequencing. We then examined the expression pattern of the protein-coding genes under these two conditions with Affymetrix Zebrafish Genome Array profiling. The potential roles of the microRNAome in the transcriptomic cold regulation in the zebrafish brain were investigated by various statistical analyses.

Results

Among the total 214 unique, mature zebrafish miRNAs deposited on the miRBase website (release 16), 175 were recovered in this study. In addition, we identified 399 novel, mature miRNAs using multiple miRNA prediction methods. We defined a set of 25 miRNAs differentially expressed under the cold and normal conditions and predicted the molecular functions and biological processes that they involve through Gene Ontology (GO) annotation of their target genes. On the other hand, microarray analysis showed that genes related to mRNA processing and response to stress were overrepresented among the up-regulated genes in cold-stress, but are not directly corresponding to any of the GO molecular functions and biological processes predicted from the differential miRNAs. Using several statistical models including a novel, network-based approach, we found that miRNAs identified in this study, either individually or together, and either directly or indirectly (i.e., mediated by transcription factors), only make minor contribution to the change in gene expression patterns under the low-temperature condition.

Conclusions

Our results suggest that the cold-stress response of mRNA expression may be governed mainly through regulatory modes other than miRNA-mediated regulation. MiRNAs in animal brains might act more as developmental regulators than thermal adaptability regulators.

Background

MicroRNAs (miRNAs) are a class of ~22 nt, non-coding RNAs present in plants and animals that play important roles in regulating gene expression by directing translational repression and/or mRNA destabilization [1-3]. Recognition and interaction between miRNA and mRNA take effect only if the 5' end seed (6-8 nt) of the miRNA is complementary to the 3' UTR site of the mRNA. A miRNA typically targets multiple, even hundreds of, genes [4], which makes miRNAs potentially global regulators responsible for spatio-temporal gene expression patterns.

Like let-7 [5,6], some miRNAs are highly conserved throughout evolution and participate in pivotal biological processes [7]. Lineage- or species-specific miRNAs may contribute to the uniqueness of organisms [8,9]. With recent advancements in miRNA microarrays and next-generation sequencing, miRNA expression profiles of many species are being generated at a rapid pace. These profiles are invaluable resources for better understanding gene expression dynamics [10-16]. Some miRNAs have been well characterized and have been shown to have important functions in normal development [17-21]. The deficient expression of some miRNAs also presents severe phenotypic abnormalities [22,23]. The role of miRNAs in responses to cellular and environmental stress has been discovered in a variety of organisms [24-28]. A few of these miRNAs have been characterized as key regulators of the stress response [25]. For example, miR-8 family miRNAs were shown to be indispensable in the response of zebrafish embryos to osmotic stress [29]. Despite this progress, most studies to date have only examined miRNAs' functions at individual level. The global contributions of microRNAomes (miRNAomes) to environmental adaptations remain largely unexplored.

The zebrafish, an important model organism, is an eurythermal fish with temperature tolerance ranging from approximately 7°C to 40°C [30], making it a useful experimental subject for exploring the molecular mechanisms underpinning thermal responses. Almost all of the previous work on the molecular underpinnings of thermal responses has focused on the fluctuated expression of protein-coding genes in response to warm or cold acclimation [31,32]. Many miRNAs have been discovered in zebrafish [33,34], and their expression profiles [35] and functions in developmental regulation have been studied [36]. Our goal was to determine how the zebrafish miRNAome responds to thermal acclimation and to assess the roles of miRNAs in thermally challenged transcriptomic plasticity. Toward this end, we systematically characterized the miRNAomes of zebrafish brains under normal and cold-acclimated conditions using Solexa sequencing. We also identified genome-wide gene expression patterns under each of these conditions using microarray profiling. We then examined the potential links between the changes in miRNA levels and the transcriptomic shifts and found that miRNAs may play minor roles in transcriptomic plasticity during cold acclimation in the zebrafish brain.

Results

Deep sequencing of the zebrafish brain small RNA and miRNA prediction

Total RNAs from brain of cold-treated (10°C, 10 days) and normal control (28°C) zebrafish were extracted and electrophoresed to enrich the 18-30 nt small RNAs and then subjected to Solexa sequencing. The numbers of high-quality reads were 9,615,499 and 14,399,653 for the normal control (NC) and cold-acclimated (WC) RNA libraries, respectively. The length distributions of the sequence reads from the two samples both peaked at 22 nt, conforming to size distribution patterns characteristic of miRNA libraries (Figure 1). After filtering short (<18 nt) and adaptor-contaminated reads, we obtained 6,473,829 (NC) and 10,799,032 (WC) redundant reads that corresponded to 349,107 (NC) and 350,474 (WC) unique reads.

thumbnailFigure 1. Length distribution of the high-quality reads in the 10-30 nt range.

We ran the SOAP2 program [37] to map the unique reads to the zebrafish genome (Danio_rerio.Zv8.56.dna.toplevel.fa.gz). Stem-loop structures were predicted from the putative miRNA precursor sequences using MIREAP (https://sourceforge.net/projects/mireap/ webcite). Overall, we identified 602 stem-loop precursors that coded for 700 miRNAs. Among them, 175 were previously known miRNAs [38], constituting over 80% of the unique miRNAs (214) deposited in the miRBase database (release 16). 525 are newly identified miRNAs (Additional file 1). To cross-validate these newly identified miRNAs, we adopted another miRNA prediction software, miRDeep [39] by invoking the quantifier.pl script (using default parameters except that -g was set to 0) in the software package. Results showed that the miRNA prediction outcomes are almost the same, independent of the prediction methods used and the miRNA expression levels are also highly consistent between the two methods (Additional files 2, 3 and 4). To further evaluate the accuracy of miRNA prediction, we examined the predicted miRNAs with a higher stringency by adding another two criteria: 1) they cannot be exon-derived (a recent study by Stark et al. [40] suggested that miRNAs lie almost exclusively in the introns and intergenic regions of fruit flies); 2) their precursors need to be recognized as genuine stem-loop structures by the MiPred [41] software. Consequently, the number of newly identified miRNAs in our dataset was reduced from 525 to 399 (See Additional file 1). However, MiPred does not ensure 100% accuracy in determining the authenticity of miRNAs. For instance, 5 of the 175 well-known zebrafish miRNAs were incorrectly assigned in this prediction. In addition, "exon-derived" miRNAs are possible to exist if the relevant exons do not code for proteins or are transcribed under alternative promoters or from alternative strand. Therefore, the number of authentic miRNAs in the zebrafish brain is very likely higher than 574 (i.e., 175+399). To verify the authenticity of the newly identified miRNAs, we selected 13 miRNAs expressed at various levels to perform RT-PCR. All of them can be successfully amplified (Additional file 5). Cloning and sequencing of the PCR products confirmed that they are in correct forms. These results showed that the newly identified miRNAs are very likely true miRNAs. Remarkably, the newly identified miRNAs are expressed at substantially lower levels than previously known miRNAs (Figure 2). This may partially explain why the new miRNAs were missed by traditional cloning and sequencing methods.

Additional file 1. Table S1. List of 525 newly identified zebrafish miRNAs.

Format: XLS Size: 155KB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Additional file 2. Figure S1. Scatter plot shows a high degree of consistency between miReap and miRDeep on miRNA identification.

Format: TIFF Size: 604KB Download fileOpen Data

Additional file 3. Data S1. Newly identified miRNAs and their precursors in normal zebrafish brain. The predicted secondary structures (bracketed notations) were produced by miRDeep2. The reads that were mapped to the precursors are aligned below, with the number of reads shown on the left.

Format: TXT Size: 538KB Download fileOpen Data

Additional file 4. Data S2. Newly identified miRNAs and their precursors in cold-acclimated zebrafish brain. The method used here is the same as stated in Additional file 3.

Format: TXT Size: 653KB Download fileOpen Data

Additional file 5. RT-PCR verification of newly identified miRNAs.

Format: DOC Size: 64KB Download file

This file can be viewed with: Microsoft Word ViewerOpen Data

thumbnailFigure 2. Violin plot showing the expression levels of novel-identified and known miRNAs. The violin plot combines box plot and kernel density trace to describe the distribution pattern of a vector of data. The four violin plots were drawn according to the expression levels of 525 newly identified and 175 known miRNAs under the normal and cold-acclimated condition, respectively. Note that if the expression level of a miRNA was 0 under certain condition, we assigned it to 1 in order to make the log-transfer of the raw data.

Solexa sequencing are known for its short reads, and errors can be pervasive in short-read sequencing [42]. To address this concern, we clustered the raw reads using the FreClu software [43], which takes the read count into consideration. The results showed that this clustering process can improve the accuracy of short reads mapping (Table 1). However, it did not improve the data quality of the predicted miRNAs (Additional file 6). Hence, we performed all of the follow-up analyses of the miRNAs identified using the original reads.

Table 1. SOAP mapping results of reads with or without clustering.

Additional file 6. Figure S2. The same pattern holds for miRNA species and expression levels with or without frequency-based reads clustering.

Format: TIFF Size: 516KB Download fileOpen Data

The miRNA expression profile in cold acclimation

Based on normalized read counts of each miRNA in the control and cold libraries, we deduced the expression changes of the brain miRNAome in response to cold acclimation. Because low-expressed sequences have more noise than high-expressed sequences, Koh et al. proposed an adaptive thresholding method using Kolomogorov-Smirnov statistics to determine the threshold value for identifying significant miRNAs [44]. We observed a similar pattern in the present study. The less abundant miRNAs showed a greater divergence than the more abundant miRNAs (Figure 3). We hypothesized that the low-expressed miRNAs showing greater divergence likely have more noise than the high-expressed ones and that the high-expressed miRNAs have a high probability of playing an important biological role. We borrowed the concept of Koh et al. to define miRNAs differentially expressed under the normal and cold conditions as the ones that have a greater than 2-fold change between the two libraries and a mean expression level greater than 64. As a result, we found 11 up-regulated and 15 down-regulated miRNAs in the cold library; a majority of them are known miRNAs (Table 2).

Table 2. Differentially expressed miRNAs of zebrafish brain under cold-acclimated and normal conditions.

thumbnailFigure 3. Scatter plot showing the miRNA expression levels of the zebrafish cold and normal brain libraries. The plot indicates much more expression noise of lowly expressed miRNAs than that of highly expressed ones. Hence, to define biologically significant miRNAs, we screened the differential genes from those with mean expression level greater than 64.

To reveal which biological processes might be regulated by the changed miRNAome, we predicted the target genes of the differential miRNAs using PITA [45]. We then projected the functions of these differential miRNAs onto the overrepresented Gene Ontology (GO) terms of their targets using the GO:TermFinder package [46] with a corrected p value of ≤0.05 as a cutoff. In total, we found 4 (dre-mir-30d-5p, dre-mir-182-5p, dre-mir-458-3p, and miR_8) and 5 miRNAs (miR_46, dre-mir-133a-2-3p, dre-mir-146a-5p, dre-mir-183-5p, and miR_8) with overrepresented biological process and molecular function GO terms, respectively (Table 3). In addition, we assembled all target genes of the up-regulated and down-regulated miRNAs and examined whether any overrepresented functions would emerge. Our results showed that genes related to "purine nucleotide binding" and "nucleotide binding" were enriched in the targets of up-regulated miRNAs.

Table 3. Differentially expressed miRNAs and their functions.

Gene expression profiling and differentially expressed genes

Gene transcription changes between the cold-acclimated and normal samples were obtained by microarray-based hybridization using three sets of RNA samples prepared from fishes of three acclimation experiments. The gene expression level was measured as the intensity of probe sets on an Affymetrix microarray. After filtering the probe sets and retrieving the map relationship between the probe sets and the genes, we obtained gene expression profile recordings of 7356 protein-coding genes. The differential genes were determined using a two-class, unpaired comparison in the SAM (version 3.0) program [47]. Using a 5% false discovery rate (FDR) cut-off, we found that 311 genes were down-regulated and 315 genes were up-regulated in the cold-acclimated samples (Additional file 7). There were 127 and 85 genes with more than 2-fold changes in mRNA abundance that were down- and up-regulated, respectively.

Additional file 7. Table S2. List of significantly up- or down-regulated genes. Genes with 2-fold changes are highlighted. Genes in red and blue are up- and down-regulated genes under the cold-acclimated condition.

Format: XLS Size: 179KB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

We then investigated the functional enrichment of these differential genes using the GO:TermFinder package [46]. As shown in Table 4, under cold stress, the genes with a functional annotation related to "mRNA processing" were enriched in the up-regulated differential genes, and genes associated with "cellular macromolecular complex subunit organization" were overrepresented in the down-regulated genes. An interesting observation here is that the identified molecular functions and biological processes changes in the cold-acclimated zebrafish brain bear no apparent overlap with those predicted from the potential targets of the differential miRNAs (see Table 3), suggesting the role of miRNAs might be minor in the transcriptomic changes in cold stress. We thus examined the possible regulation roles of the microRNAome in cold adaptation by other statistical approaches.

Table 4. Functional enrichment of differentially expressed mRNAs.

Verification of the expression pattern of the protein-coding genes and miRNAs by quantitative RT-PCR (q RT-PCR)

The reliability of statistical analyses depends on the quality of the microarray and solexa profiling data. To evaluate the accuracy of the microarray screening and the miRNA profiling results, genes and miRNAs with change direction either down-regulated (plp1b, h3f3c and dre-mir-24-4-3p) or up-regulated (emx2, ikzf5, dazap, mmp13, zgc:112425, taldo1, and dre-mir-152-3p) were selected for quantitative RT-PCR verification. Q-PCR showed that the direction of expressional change for each selected gene correlated perfectly with the results of the microarray screening or solexa sequencing (Figure 4). Although the changed folds of expression in a few genes (mmp13a, plp1b) and the dre-mir-152-3p showed some discrepancy between the q-PCR results and the '-omic' (microarray and solexa sequencing) profiling data, the other genes/miRNAs are in good agreement. Since we evaluate only the direction of change of the genes and miRNAs when delineating the potential relationship between miRNA-mediated regulation and transcriptomic shifts (see later sections), the gene and miRNA expression profiles obtained from the '-omic' profiling approaches are generally useful for statistical analyses.

thumbnailFigure 4. Fold change of expression level for selected genes and miRNAs measured by q RT-PCR. Quantitative RT-PCR shows that the direction of expressional change under the cold condition of the selected genes and miRNAs correlated nicely with the results obtained from the "-omics" data (Additional file 7 and Table 2).

Statistical tests for the potential correlation between the expression levels of miRNAs and protein-coding genes

Mounting evidence has shown that miRNAs can orchestrate the expression of mRNAs and that the well-balanced interplay between the miRNAome and protein-coding gene expression plays a crucial role in an organism's development [48-50]. To investigate the influence of miRNAs on the global expression of protein-coding genes in cold stress, we adopted four different approaches. First, we estimated the effect of miRNAs on differentially regulated genes using Fisher's exact test. Using PITA target prediction methods, we obtained 18665 pairs of miRNA-target interactions in total. Our results showed only one miRNA (dre-mir-34-3p) with target genes overrepresented in the down-regulated gene set. No miRNAs showed targets enrichment or depletion in the up-regulated gene set (FDR < 5%). The miRNA dre-mir-34-3p was expressed at very low levels (29 in the cold library and 20 in the control library) in the brain and was not significantly different between the two conditions. We speculate that it is not probable that dre-mir-34-3p plays an important role in the changed gene expression patterns observed under cold stress. To address concerns over potential target prediction bias by PITA, we performed target prediction by miRanda software [51]. Under this methodology, we obtained 56670 pairs of miRNA-target interactions, which is about 3-fold of the number (18665) obtained by PITA prediction method. We then examined the composition of targets predicted by the two methods for each miRNA. Our results showed that 615 miRNAs shared at least one predicted target between PITA and miRanda; however, as a whole, the overlap is relatively low (Additional file 8, Table S3-S6). For example, among a group of 370 miRNAs for which each has more than ten targets by whatever target prediction methods, the median value of the proportion of overlap is only 0.396 (Additional file 9, Figure S3). In average, each miRNA targets ~80 protein-coding genes and each gene is likely to be regulated by 7-8 miRNAs in return. Based on these interaction data, our statistical analyses recaptured the fore-mentioned miRNA-mRNA interacting pattern, i.e., no miRNAs targeted over- or under-represented genes in either up- or down-regulated gene set.

Additional file 8. Table S3. List of miRNA-target pairs predicted by PITA. Table S4. List of miRNA-target pairs predicted by miRanda. Table S5. List of miRNA-target pairs predicted by PITA and miRanda. Table S6. Summary of target overlap between PITA and miRanda.

Format: XLS Size: 4.2MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Additional file 9. Figure S3. Histogram showing the distribution of miRNA targets commonly predicted by PITA and miRanda. 370 miRNAs were counted, each of which targeted more than 10 genes by whatever prediction methods used. The degree of overlap is calculated as the ratio of the number of shared targets to the minimal value of the number of targets predicted by PITA and miRanda, separately.

Format: TIFF Size: 126KB Download fileOpen Data

Recently, Chen et al. [52] proposed a novel concept for the regulatory effect score (RE score) to measure the inhibition capability of a miRNA on its targets. In principle, the RE score is calculated as the average difference in expression levels between a target and a non-target group. By calculating the RE scores on the PITA and miRanda methodology, respectively, we checked whether some of the miRNAs would changed regulatory effects on their targets between the two conditions. Consistent with previous results, we again found that no miRNAs exhibited a significant change in their regulatory effect after correcting the raw t-test p-values by the multiple hypothesis test (FDR < 0.05) whatever the target prediction method was used. Taken together, our results indicate that individual miRNAs are unlikely to be major factors leading to the changed transcriptomic phenotype of zebrafish brain under cold stress.

To further reduce the concerns over target-prediction bias induced errors, we applied the Sylamer algorithm [53] to re-examine the potential relationship between the miRNAome and transcriptome under cold acclimation. We sorted the protein-coding genes from the most significantly up-regulated to the most down-regulated. We then calculated the distribution of the enrichment/depletion of words in the 3' UTR sequences that could be complementary to the seed sequences of miRNAs. Since positions 1-7 and 2-8 of the 5' end of a mature miRNA are the most important sequences for target recognition [54], we took the 7 nt region as a miRNA seed in this study and examined whether some over- or under-represented words corresponded to certain miRNAs.

As shown in Figure 5, there are no obvious peaks in the region containing the most down-regulated genes, suggesting that the analysis failed to identify any miRNAs being the major and direct factors contribute to the decreased expression levels of the down-regulated gene sets; on the other hand, few than ten words seem to be overrepresented in the 3' UTR regions of the up-regulated genes. The three most overrepresented words are TAATAGA (P = 3.86 × 10-10), GCGAATT (P = 1.58 × 10-9) and GTTCACG (P = 4.58 × 10-9). However, these words are not complementary to any seed region of the profiled zebrafish miRNAs. We then examined the composition of the top 10 most enriched words and found only one, GCTAGTT (P = 5.53 × 10-8), that was complementary to the seed region of the newly identified miR_381 (AACTAGCAGCTGTTGACATCCA). However, this miRNA was hardly expressed in the brain (1 in the control library and 4 in the cold library), suggesting that it was unlikely to be the major factor responsible for the changed gene expression patterns observed.

thumbnailFigure 5. The word enrichment and depletion pattern in mRNA 3' UTR. Genes are sorted by their expression divergence between normal and cold-acclimated zebrafish (from the most up-regulated to the most down-regulated). Probability values with a plus-sign are enriched; those with a minus-sign are depleted.

To systematically investigate the extent to which miRNAome-mediated regulation could have influenced the transcriptome under cold-acclimation, we designed a novel, network-based methodology to parse the subtle regulatory effects of a group of miRNAs on a group of protein-coding genes at the system level. Based on the predicted miRNA-target relationship predicted from the PITA method, we first built a miRNA-mRNA network of 32 differential miRNAs and 626 mRNAs (311 down-regulated and 325 up-regulated). The miRNA-target relationship showed that the majority of the differential gene transcripts do not possess any potential targeting sites in the 3' UTRs that correspond to the differentiated miRNAs. In addition, approximately half of the differential miRNAs do not target any genes in the differentiated gene list. The resultant miRNA-mRNA network involved only 97 nodes (82 mRNAs and 15 miRNAs) and 97 edges when the isolated nodes were separated. As seen in Figure 6, the miRNA-mRNA interaction networks are composed of six unconnected components. One component contains 60 of the 97 nodes, whereas the others are formed by only 1 or 2 miRNAs situated in the network center. The node with the highest number of edges is the up-regulated miRNA miR_46 that was predicted to target 18 (11 up-regulated and 7 down-regulated) gene transcripts. Among the 82 miRNA-targeted differentially expressed genes, over four-fifths (67/82) are predicted to be targeted by only 1 miRNA and one-fifth by 2 or more. We then computed the "regulatory density" (RD) of the network by formula (1) (see Methods). We developed this formula and define it as a proxy for the regulatory potential of a group of miRNAs upon gene transcripts by calculating the proportion of coherent over incoherent regulatory links within a miRNA-mRNA network. Our analysis showed that the observed RD was 0.1134. However, the values of the 95th- and 99th-quantiles were 0.1546 and 0.1959, respectively, for the 1000 cohorts of random networks that were constructed by shuffling the original miRNA-mRNA network while maintaining the degree distribution. We further analyzed the topologic property of a similar network derived by using the miRanda prediction data. The produced miRNA-mRNA network consisted of 203 edges that link 25 and 161 differential miRNAs and protein-coding genes, respectively. As seen in Additional file 10, a majority of the mRNAs are targeted by a few miRNAs; and only a few miRNAs have >5 target genes. Calculation showed that the RD score of this network is 0.0148, again not significantly greater than random expectation (empirical P-value: 0.269). Taken together, results from all the above statistical analyses suggested that the miRNAome is unlikely to be the major regulatory component responsible for the changed expression profiles of the protein-coding genes under cold stress.

Additional file 10. Figure S4. Graphic illustration of the regulatory relationship between the differential miRNAs and mRNAs as predicted by miRanda. The style of this graph is the same as Figure 6.

Format: TIFF Size: 1.1MB Download fileOpen Data

thumbnailFigure 6. Graphic illustration of the regulatory relationship between the differential miRNAs and mRNAs predicted by PITA. Triangles, circles and diamonds denote miRNAs, non-TF genes and TF genes, respectively. The nodes in red and blue denote up- and down-regulated miRNAs or mRNAs, respectively. The edges in dark and yellow denote coherent and incoherent regulation, respectively. The size of a node is proportional to the number of its partners. Of nodes with degree > 5, the associated name is shown.

Mathematical models [55] and statistical associations [56,57] have shown that miRNAs frequently intertwine with transcription factors (TFs) in regulatory circuits and play essential roles in signal transduction. Recently, Tu et al. [58] proposed a two-layer network wherein miRNA-initiated regulation was mediated and propagated through the miRNAs' targeted transcription factors. We investigated whether the miRNAs had indirectly, through their transcription factors, led to the observed transcriptomic plasticity. We first filtered the transcription factors from the arrayed gene list according to the GO annotation. We treated the genes with any one of the GO terms, "DNA binding", "transcription factor activity", "transcription activator activity" or "transcription repressor activity", as putative transcription factors. Accordingly, we predicted 796 TF-coding genes among the 7356 protein-coding genes. Notably, only 46 of these TFs showed differential expression in cold acclimation, and they were significantly underrepresented in the differential expressed gene list (46 out of 626 differential genes versus 796 TFs in the 7356 arrayed genes, Fisher's exact test: P < 0.01). Among the differential TF genes, only 6 were predicted to be the targets of 5 differential miRNAs (the ◆-▲ pair in Figure 6). Of the 6 TFs, 3 showed coherent regulatory direction with their respective miRNA partners (miR_13→emx2; dre-mir-183-5p→ikzf5; miR_46, miR_48→h3f3c. see Figure 6). Given that the gene emx2 represented slightly up-regulatoin (fold-change: 1.214; q-value: 4.376%) under cold-acclimated condition, we speculate that it was not a driven gene responsible for the changed mRNA expression pattern. There is no evidence that the ikzf5 (fold-change: 1.994; q-value: 0) and h3f3c (fold-change: 2.245; q-value: 2.144%) are related to stress response when we examined their annotations in multiple databases and in the literature. Additionally, almost no differential TF-coding genes were targeted by differential miRNAs. Taken together, we conclude that the mRNA phenotypic plasticity seen under the cold environment was not dominantly modulated by miRNAs, even if these factors were directly or indirectly mediated by transcriptional factors.

Discussion

The repertoires of miRNAs have been constantly revised and updated with the advance of in silico prediction methods and more intensive cloning and sequencing efforts in broader tissues and developmental stages. The number of unique, mature human miRNAs is now up to 1223 in the miRBase database (version 16). In contrast, there are only 214 mature miRNAs now registered for zebrafish. This relatively low number could be the result of genomic between-species differences or incomplete miRNA identification in zebrafish. By Solexa sequencing, we comprehensively profiled the miRNA repertoire in the zebrafish brain. The number of predicted miRNAs in this tissue reached 700 (574 in high-stringency criteria) in this study and greatly expanded the miRNAome in zebrafish. Of note is that the expression levels of the newly identified miRNAs in the zebrafish brain are mostly low. The authenticity of these newly identified miRNAs remains to be further investigated to determine their biological functions. However, a population analysis on 7 randomly selected novel miRNAs currently carried out in the laboratory showed that almost half of them are subjected to purified selection suggesting that many newly identified miRNAs may have high potentials to be functional in vivo despite their low expression levels. It is possible that some of these newly identified, low-expressed miRNAs could be expressed at higher levels in tissues other than the brain and play roles in gene regulation. The miRNAome of the zebrafish is much larger than those of the lower chordates such as the amphioxus species, which had been determined by small RNA library sequencing [59] and solexa sequencing [13]. This suggests the expansion of the miRNA regulome in the evolution of vertebrates.

To study the function of the miRNAome in environmental adaptation, we designed the cold-acclimation experiments in this study. We profiled both the miRNAome and the transcriptome of the brain under low- and normal-temperature conditions. We adopted four approaches to examine how miRNAs, either individually or globally, could have influenced the transcriptomic plasticity under cold conditions. Using statistical methods we first examined the potential effects of individual miRNAs on gene expression patterns. We found one or two miRNAs exhibiting regulation signatures on mRNA expression based on enrichment/depletion analysis of target genes among the differential gene set. However, we argue that these miRNAs are unlikely to have been responsible for the different expression patterns after further investigation of their expression levels. Our results show that at the individual level, the miRNAs were not a main factor shaping the mRNA expression patterns observed under the cold-acclimated conditions. We then attempted to clarify whether the cold-stress response of mRNA expression is governed by miRNAs indirectly and is mediated by TFs. TF-mediated transcriptional regulation is thought to be stronger than the post-transcriptional regulation usually mediated by miRNAs [60]. For this purpose, we assigned a set of TFs according the GO annotation and found that TFs were significantly underrepresented in the differential protein-coding genes; there were only 6 TFs targeted by differential miRNAs. Again, our results suggest that there is hardly any influence from miRNAs on the changed mRNA expression patterns. A recent study by Nunez-lglesias et al. [50] suggested that a subtle regulation relationship between miRNAs and their targets can be observed in the global miRNAs-mRNA network using more sensitive experimental designs and a systems biology approach. Aiming to discern the subtle regulatory effects, we employed a novel network-based analysis. We found that the miRNA-mRNA network (formed by the regulatory relationship between differential miRNAs and differential mRNAs) does not possess a significantly higher percentage of coherent edges than the random expectation. Furthermore, this regulatory network consists of only 3 pairs of coherent edges between miRNAs and TF genes. Therefore, our data also rule out the possibility that miRNAs as major factors indirectly modulate mRNA expression in the brain of cold-acclimated zebrafish. One might argue that many TFs are often expressed in low abundance and that some of them might go beyond the detective sensitivity of the current Affymetrix array screening methodology, causing the observed underrepresentation of the TFs in the differential gene set. This situation might eventually lead to incomplete evaluation of the regulatory effects of miRNAs on the transcriptome. To evaluate this possibility, we calculated the mean expression level of the 796 predicted TFs and the non-TF genes in the array. As seen from figure 7, though the TFs showed slightly lower average expression levels than the non-TF genes, the distribution pattern of expression intensity for TFs as compared with non-TFs within the differential genes is very similar as that within the whole set of the arrayed genes, indicating that the relatively low expression of TFs did not hinder the screening of the differential TF genes in this study. Evidence from various developmental biology studies has been continuing emerged favoring the hypothesis that miRNAs are capable of imparting robustness on gene regulation networks that involve several classes of recurrent circuits composed of miRNAs and their corresponding targets [61]. For example, a study by Hilgers et al. [62] revealed that miR-263a/b plays an important role in protecting Drosophila mechanosensory bristles from apoptosis and ensures the developmental robustness of this organ. In our study, we found miRNAs are minor factors in transcriptome changes in the cold acclimated zebrafish brain. This difference may reflect that miRNAs in animal brains act more as developmental regulators than adaptability regulators.

thumbnailFigure 7. The expression pattern of TFs and non-TFs within differential or non-differential gene sets. The boxplot shows that the TFs expressed at a slightly lower level than the non-TFs. However, the expression pattern for TFs as compared to non-TFs within the differential genes is very similar as that within the whole set of the arrayed genes, indicating that the relatively low expression of TFs did not influence the screening of the differential TF genes in this study. DEG: differential genes; non-DEG: non-differential genes.

Thermally invoked responses can be generally classified as cellular stress responses, homeostatic responses or body remodeling-related responses that correspond to short-term, mid-term or long-term acclimation processes, respectively [63-65]. Our samples were taken after 10-day cold acclimation, which can be classified as mid-term acclimation. We reason that the current results reflect only the gene regulatory patterns of the zebrafish brain at this temporal scale of acclimation. The involvement of miRNAs as a pivotal regulator in an earlier or later cold adaptation phase remains to be further investigated.

Persistent efforts have identified and characterized some cold stress-induced protein-coding genes [66,67]; among them, rbm3, which encodes a glycine-rich RNA binding protein, has been suggested to regulate global protein synthesis by affecting miRNA expression levels [68]. In the present study, we also observed two homologous genes in zebrafish, cirbp (fold-change: 2.244; q-value: 0) and zgc:112425 (fold-change: 1.961; q-value: 1.598%), with significantly up-regulated expression levels under cold-acclimated condition.

Lopez-Maury et al. [69] have pointed out that post-transcriptional regulatory mechanisms, compared to transcriptional mechanisms, are largely unknown under stress conditions. To our knowledge, there are few studies jointly analyzing miRNAs and mRNA profiles to explore the roles of miRNAs on the cold adaptation of fish. Such a study would extend our understanding of miRNAs and environmental adaptation and undoubtedly lead to a more challenging and interesting question: whether and how miRNAs have integrated into cold-adaptive processes and contributed to the unique physiology of fishes endemic to extreme temperature conditions.

Conclusions

In this study, we systematically characterized the miRNAome of zebrafish brain under normal and cold-acclimated conditions using Solexa sequencing. Our data has greatly extended the zebrafish miRNAs repertory from the 214 unique, mature miRNAs (miRBase databse, release 16 version) to more than 574 ones. Besides, our results showed that the miRNAs, either individually or together, either directly or indirectly play only a minor role in transcriptomic plasticity during cold acclimation in the zebrafish brain, suggesting that the cold-stress response of mRNA expression may be governed mainly through regulatory modes other than post-transcriptional regulation by miRNAs. MiRNAs in animal brains may act more as developmental regulators than thermal adaptability regulators.

Methods

Zebrafish treatment and brain RNA preparation

Adult zebrafish (Danio rerio) were alternatively divided into two groups: a control group and an experimental group. In each group, the sex ratio was 1:1. The control group was raised under the normal laboratory breeding temperature (28°C), and the experimental group was subjected to cold (10°C) stress. An abrupt transfer of zebrafish from normal to cold (10°C) environments will cause most of them to die. To prevent this, the temperature of the cold group was decreased from 28°C to 10°C at a rate of one degree per hour. After 10 days of cold adaption, the brains of the zebrafish were dissected, and total RNAs were extracted from these brains using the total RNA isolation kit (Promega, Madison, USA). All experimental research on animals was conducted with the approval of the Animal Research Ethics Committee of our institute.

Microarray profile

To evaluate the genome-wide change of gene transcription in cold stress, we performed hybridization analysis using the Zebrafish Genome Array (Affymetrix) as previously described [70]. Briefly, mRNAs were isolated from 10 μg of total RNAs derived from pooled zebrafish brains of acclimated and control samples using the Eukaryotic Poly-A RNA Control Kit (Affymetrix) in an RNase-free environment. Subsequently, mRNAs were subjected to double-strand cDNA synthesis using a one-cycle cDNA synthesis kit (Affymetrix). The resulting double-strand cDNAs were in vitro transcribed into biotin-labeling cRNAs using the IVT Labeling kit (Affymetrix). After that, biotin-labeling cRNAs were fragmented by heating and hybridized with the Affymetrix Zebrafish Genome Arrays. Finally, the signal intensity of the chip was scanned using the Gene Chip Scanner 3000 (Affymetrix) and analyzed with the Gene Chip Operating Software (Affymetrix). The data derived from the genome-wide expressional analysis was deposited in the GEO bank with the accession number GSE32092.

Small-RNA deep sequencing

For small-RNA Solexa sequencing, 20 μg of total RNAs were subjected to electrophoresis on a polyacrylamide gel under denaturing conditions. The small RNAs ranging from 18-30 nt were excised and purified. The resulting small RNAs were sequentially ligated to a 5' adaptor and a 3' adaptor. The products were then reverse transcribed (RT) and amplified by polymerase chain reaction (PCR). The RT-PCR products of approximately 90 bp (small RNA + adaptors) were isolated and subjected to sequencing analysis using the Illumina Genome Analyzer (Illumina, San Diego, USA) according to the manufacturer's instructions.

MiRNA identification and annotation

High-quality, small-RNA reads larger than 18 nt were extracted from raw reads and mapped to zebrafish genome sequences (Danio_rerio.Zv8.56.dna.toplevel.fa.gz) using SOAP2 [37]. The miRNAs and their precursor structures were predicted by MIREAP (https://sourceforge.net/projects/mireap/ webcite) under the default settings using the perfect-matched reads. Like the miRDeep software [39], MIREAP integrates miRNA biogenesis, sequencing depth and structural features to identify miRNAs and further interrogate their expression level from deep sequenced small RNA libraries. Stem-loop hairpins were retained only when they comply with: 1) the mature miRNAs-associated reads are mapped in the arm region of the precursors; 2) the free energy of the secondary structure calculated by RNAfold [71] is lower than -18 kcal/mol. In detail, the potential mature miRNA is defined as the most abundant read sequence that aligns to the potential precursor sequence; the expression level of a miRNA is then specified as the sum of an ensemble of reads that align with the potential mature molecules, allowing three nucleotides sliding beyond the position of the potential mature miRNA at the 5' end. We then performed blast analyses on the predicted miRNAs precursors and mature miRNAs against sequences deposited in the miRBase (release 16) database to annotate these miRNAs. Novel miRNAs were further tagged as "exon-derived", "intron-derived" or "intergenic" miRNAs according to the loci of their precursors on the genome. In addition, MiPred [41] was used to distinguish between genuine and pseudo miRNA precursors. This software incorporates the local, contiguous structure-sequence composition and the minimum free energy (MFE) of the secondary structure to eliminate potential false positive predictions.

Frequency-based small-RNA reads clustering

Following the method proposed by Qu et al. [43], the raw reads were clustered according to their read numbers. After clustering, the accuracy of mapping short reads can be improved.

Gene expression validation by q RT-PCR analysis

Total RNA (8 μg) from the cold and normal fishes was respectively reverse transcribed by Takara PrimeScript RT-PCR Kit (Takara, CN) with OligodT primers, according to the protocol the manufacturer recommended. Another 4 μg total RNA from the two treatments was used for small RNA (< 200 bp) isolation by Ambion mirVana miRNA Isolation Kit. Poly (A) tail was added to the small RNAs by Ambion Poly (A) Polymerase (Ambion, CN). 2 pairs of primers for each selected gene or miRNA were designed and tested for amplification efficacy and specificity by agarose electrophoretic analysis of the PCR products. The primer pair with the best specificity was selected for further use in q-PCR. Q-PCR was run by using the SYBR RT-PCR kit in a Bio-rad CFX 96 (Bio-rad) machine and results were analyzed by the CFX Manager software. The PCR thermal cycling was performed using 45 cycles of 95°C for 15 seconds, 72°C for 30 seconds and 60°C for 30 seconds. The PCR reaction for each gene was performed in triplet with either the housekeeping gene Beta-actin (for arrayed genes) or U6 RNA (for miRNAs) as control.

Word enrichment/depletion in 3' UTR

The mRNAs were sorted according to their significance of expression levels and their corresponding 3' UTRs. The longest sequences were retrieved from the Ensembl website using Biomart tools [49]. The word enrichment/depletion patterns of these genes were generated by Sylamer [53].

Regulatory effect score

The regulatory effect score (RE score) was defined as the inhibitory effect of a miRNA on the mRNAs in a sample, according to the average difference in expression of its target versus non-target mRNAs [52]. For each miRNA, we first calculated the rank-based RE score for every sample and then compared the difference between the normal and cold-acclimated groups using a two-sample t-test. The FDR values were calculated using the R multtest [72] package.

GO enrichment test

For the list of mRNAs, we tested whether each had enriched GO terms in biological processes and molecular functions using the GO:TermFinder package [46]. The GO annotation file was downloaded from http://zfin.org/zf_info/downloads.html#go/on webcite January 29, 2010. The GO ontology file was downloaded from http://www.geneontology.org/ webcite on January 23, 2008.

miRNA target predictions

Currently, there are several popular tools for predicting miRNA targets, including PITA [45], PicTar [73], miRanda [51], mirWIP [74] and TargetScan [75]. These tools can be generally divided into two classes, methods that integrate the information of sequence conservation across multiple genomes and those that do not take this into account. Because of the limited genome data for zebrafish relatives, we predicted the miRNA targets by two methods, PITA and miRanda, independently. PITA scores the miRNA-target interactions based on a thermodynamic model. A salient feature of this program is that it considers the site accessibility of the surrounding mRNA sequence but disregards sequence conservation between genomes. Thus, non-conserved or even species-specific target genes could be identified using this method. For each miRNA, the target genes are defined as those predicted to contain at least one binding site (△△G ≤ -10) at their 3' UTR. We also ran miRanda (August 2010 Release) to predict the target genes with default parameters except that the energy parameter is set to -20.

Network analysis

The miRNA-mRNA relationship is modeled as a bipartite where an edge links a miRNA and an mRNA with a predicted binding site in its 3' UTR corresponding to the miRNA. We defined coherent and incoherent edges as those for which their associated miRNAs and mRNAs represented the opposite and parallel expression change under cold stress, respectively. We then calculated a regulatory density (RD) for a group of miRNAs and mRNAs within a network according to the following formula:

<a onClick="popup('http://www.biomedcentral.com/1471-2164/12/605/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2164/12/605/mathml/M1">View MathML</a>

(1)

where #CE (#IE) denotes the number of coherent (incoherent) edges.

To estimate the significance of the observed RE, we generated 1,000 cohorts of random networks and recalculated the RE by shuffling the original network while keeping the degree preservation according to the published method [76].

Authors' contributions

RY, ZD and LC conceived of and designed the study; ZD and SC performed the experiments; RY analyzed the data; RY and LC wrote the paper. All authors read and approved the final manuscript.

Acknowledgements

We are indebted to two anonymous reviewers for their constructive comments and suggestions that have substantially improved the paper. This work is financially supported by grants: 2010CB126304 from Ministry of Science and Technology of China, and 30910103906, 30971625 from Natural Science Foundation of China to Liangbiao Chen.

References

  1. Lim LP, Lau NC, Garrett-Engele P, Grimson A, Schelter JM, Castle J, Bartel DP, Linsley PS, Johnson JM: Microarray analysis shows that some microRNAs downregulate large numbers of target mRNAs.

    Nature 2005, 433(7027):769-773. OpenURL

  2. Baek D, Villen J, Shin C, Camargo FD, Gygi SP, Bartel DP: The impact of microRNAs on protein output.

    Nature 2008, 455(7209):64-71. OpenURL

  3. Hendrickson DG, Hogan DJ, McCullough HL, Myers JW, Herschlag D, Ferrell JE, Brown PO: Concordant regulation of translation and mRNA abundance for hundreds of targets of a human microRNA.

    PLoS Biol 2009, 7(11):e1000238. OpenURL

  4. Bartel DP: MicroRNAs: target recognition and regulatory functions.

    Cell 2009, 136(2):215-233. OpenURL

  5. Lee RC, Feinbaum RL, Ambros V: The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14.

    Cell 1993, 75(5):843-854. OpenURL

  6. Roush S, Slack FJ: The let-7 family of microRNAs.

    Trends Cell Biol 2008, 18(10):505-516. OpenURL

  7. Cordes KR, Sheehy NT, White MP, Berry EC, Morton SU, Muth AN, Lee TH, Miano JM, Ivey KN, Srivastava D: miR-145 and miR-143 regulate smooth muscle cell fate and plasticity.

    Nature 2009, 460(7256):705-710. OpenURL

  8. Zhang R, Wang YQ, Su B: Molecular evolution of a primate-specific microRNA family.

    Mol Biol Evol 2008, 25(7):1493-1502. OpenURL

  9. Li J, Liu Y, Dong D, Zhang Z: Evolution of an X-linked primate-specific microRNA cluster.

    Mol Biol Evol 2009. OpenURL

  10. Landgraf P, Rusu M, Sheridan R, Sewer A, Iovino N, Aravin A, Pfeffer S, Rice A, Kamphorst AO, Landthaler M, et al.: A mammalian microRNA expression atlas based on small RNA library sequencing.

    Cell 2007, 129(7):1401-1414. OpenURL

  11. Liang Y, Ridzon D, Wong L, Chen C: Characterization of microRNA expression profiles in normal human tissues.

    BMC Genomics 2007, 8:166. OpenURL

  12. Hua YJ, Tang ZY, Tu K, Zhu L, Li YX, Xie L, Xiao HS: Identification and target prediction of miRNAs specifically expressed in rat neural tissue.

    BMC Genomics 2009, 10:214. OpenURL

  13. Chen X, Li Q, Wang J, Guo X, Jiang X, Ren Z, Weng C, Sun G, Wang X, Liu Y, et al.: Identification and characterization of novel amphioxus microRNAs by Solexa sequencing.

    Genome Biol 2009, 10(7):R78. OpenURL

  14. Chen PY, Manninga H, Slanchev K, Chien M, Russo JJ, Ju J, Sheridan R, John B, Marks DS, Gaidatzis D, et al.: The developmental miRNA profiles of zebrafish as determined by small RNA cloning.

    Genes Dev 2005, 19(11):1288-1293. OpenURL

  15. Zhang J, Xu Y, Huan Q, Chong K: Deep sequencing of Brachypodium small RNAs at the global genome level identifies microRNAs involved in cold stress response.

    BMC Genomics 2009, 10:449. OpenURL

  16. Jagadeeswaran G, Zheng Y, Sumathipala N, Jiang H, Arese E, Soulages JL, Zhang W, Sunkar R: Deep sequencing of small RNA libraries reveals dynamic regulation of conserved and novel microRNAs and microRNA-stars during silkworm development.

    BMC Genomics 11(1):52. OpenURL

  17. Yoo AS, Staahl BT, Chen L, Crabtree GR: MicroRNA-mediated switching of chromatin-remodelling complexes in neural development.

    Nature 2009, 460(7255):642-646. OpenURL

  18. Johnston RJ Jr, Chang S, Etchberger JF, Ortiz CO, Hobert O: MicroRNAs acting in a double-negative feedback loop to control a neuronal cell fate decision.

    Proc Natl Acad Sci USA 2005, 102(35):12449-12454. OpenURL

  19. Schratt GM, Tuebing F, Nigh EA, Kane CG, Sabatini ME, Kiebler M, Greenberg ME: A brain-specific microRNA regulates dendritic spine development.

    Nature 2006, 439(7074):283-289. OpenURL

  20. Conaco C, Otto S, Han JJ, Mandel G: Reciprocal actions of REST and a microRNA promote neuronal identity.

    Proc Natl Acad Sci USA 2006, 103(7):2422-2427. OpenURL

  21. Makeyev EV, Zhang J, Carrasco MA, Maniatis T: The MicroRNA miR-124 promotes neuronal differentiation by triggering brain-specific alternative pre-mRNA splicing.

    Mol Cell 2007, 27(3):435-448. OpenURL

  22. Giraldez AJ, Cinalli RM, Glasner ME, Enright AJ, Thomson JM, Baskerville S, Hammond SM, Bartel DP, Schier AF: MicroRNAs regulate brain morphogenesis in zebrafish.

    Science 2005, 308(5723):833-838. OpenURL

  23. Hansen T, Olsen L, Lindow M, Jakobsen KD, Ullum H, Jonsson E, Andreassen OA, Djurovic S, Melle I, Agartz I, et al.: Brain expressed microRNAs implicated in schizophrenia etiology.

    PLoS One 2007, 2(9):e873. OpenURL

  24. Jian X, Zhang L, Li G, Wang X, Cao X, Fang X, Chen F: Identification of novel stress-regulated microRNAs from Oryza sativa L.

    Genomics 2010, 95(1):47-55. OpenURL

  25. Zhao M, Ding H, Zhu JK, Zhang F, Li WX: Involvement of miR169 in the nitrogen-starvation responses in Arabidopsis.

    New Phytol 2011. OpenURL

  26. Frazier TP, Sun G, Burklew CE, Zhang B: Salt and Drought Stresses Induce the Aberrant Expression of microRNA Genes in Tobacco.

    Mol Biotechnol 2011. OpenURL

  27. Biggar KK, Storey KB: The emerging roles of microRNAs in the molecular responses of metabolic rate depression.

    J Mol Cell Biol 2010. OpenURL

  28. Simone NL, Soule BP, Ly D, Saleh AD, Savage JE, Degraff W, Cook J, Harris CC, Gius D, Mitchell JB: Ionizing radiation-induced oxidative stress alters miRNA expression.

    PLoS One 2009, 4(7):e6377. OpenURL

  29. Flynt AS, Thatcher EJ, Burkewitz K, Li N, Liu Y, Patton JG: miR-8 microRNAs regulate the response to osmotic stress in zebrafish embryos.

    J Cell Biol 2009, 185(1):115-127. OpenURL

  30. Cortemeglia , Cheryl , Beitinger , Thomas LG: Temperature tolerances of wild-type and red transgenic zebra danios. Volume 134. Bethesda, MD, ETATS-UNIS: American Fisheries Society; 2005::7.

  31. Vergauwen L, Benoot D, Blust R, Knapen D: Long-term warm or cold acclimation elicits a specific transcriptional response and affects energy metabolism in zebrafish.

    Comparative Biochemistry and Physiology - Part A: Molecular & Integrative Physiology 2010, 157(2):149-157. OpenURL

  32. Malek RL, Sajadi H, Abraham J, Grundy MA, Gerhard GS: The effects of temperature reduction on gene expression and oxidative stress in skeletal muscle from adult zebrafish.

    Comp Biochem Physiol C Toxicol Pharmacol 2004, 138(3):363-373. OpenURL

  33. Soares AR, Pereira PM, Santos B, Egas C, Gomes AC, Arrais J, Oliveira JL, Moura GR, Santos MA: Parallel DNA pyrosequencing unveils new zebrafish microRNAs.

    BMC Genomics 2009, 10:195. OpenURL

  34. Thatcher EJ, Bond J, Paydar I, Patton JG: Genomic organization of zebrafish microRNAs.

    BMC Genomics 2008, 9:253. OpenURL

  35. Wienholds E, Kloosterman WP, Miska E, Alvarez-Saavedra E, Berezikov E, de Bruijn E, Horvitz HR, Kauppinen S, Plasterk RH: MicroRNA expression in zebrafish embryonic development.

    Science 2005, 309(5732):310-311. OpenURL

  36. Flynt AS, Li N, Thatcher EJ, Solnica-Krezel L, Patton JG: Zebrafish miR-214 modulates Hedgehog signaling to specify muscle cell fate.

    Nat Genet 2007, 39(2):259-263. OpenURL

  37. Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, Wang J: SOAP2: an improved ultrafast tool for short read alignment.

    Bioinformatics 2009, 25(15):1966-1967. OpenURL

  38. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ: miRBase: tools for microRNA genomics.

    Nucleic Acids Res 2008, (36 Database):D154-158. OpenURL

  39. Friedlander MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, Rajewsky N: Discovering microRNAs from deep sequencing data using miRDeep.

    Nat Biotechnol 2008, 26(4):407-415. OpenURL

  40. Stark A, Kheradpour P, Parts L, Brennecke J, Hodges E, Hannon GJ, Kellis M: Systematic discovery and characterization of fly microRNAs using 12 Drosophila genomes.

    Genome Res 2007, 17(12):1865-1879. OpenURL

  41. Jiang P, Wu H, Wang W, Ma W, Sun X, Lu Z: MiPred: classification of real and pseudo microRNA precursors using random forest prediction model with combined features.

    Nucleic Acids Res 2007, (35 Web Server):W339-344. OpenURL

  42. Erlich Y, Mitra PP, delaBastide M, McCombie WR, Hannon GJ: Alta-Cyclic: a self-optimizing base caller for next-generation sequencing.

    Nat Methods 2008, 5(8):679-682. OpenURL

  43. Qu W, Hashimoto S, Morishita S: Efficient frequency-based de novo short-read clustering for error trimming in next-generation sequencing.

    Genome Res 2009, 19(7):1309-1315. OpenURL

  44. Koh W, Sheng CT, Tan B, Lee QY, Kuznetsov V, Kiang LS, Tanavde V: Analysis of deep sequencing microRNA expression profile from human embryonic stem cells derived mesenchymal stem cells reveals possible role of let-7 microRNA family in downstream targeting of hepatic nuclear factor 4 alpha.

    BMC Genomics 11(Suppl 1):S6. OpenURL

  45. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E: The role of site accessibility in microRNA target recognition.

    Nat Genet 2007, 39(10):1278-1284. OpenURL

  46. Boyle EI, Weng S, Gollub J, Jin H, Botstein D, Cherry JM, Sherlock G: GO::TermFinder--open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes.

    Bioinformatics 2004, 20(18):3710-3715. OpenURL

  47. Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response.

    Proc Natl Acad Sci USA 2001, 98(9):5116-5121. OpenURL

  48. Nielsen JA, Lau P, Maric D, Barker JL, Hudson LD: Integrating microRNA and mRNA expression profiles of neuronal progenitors to identify regulatory networks underlying the onset of cortical neurogenesis.

    BMC Neurosci 2009, 10:98. OpenURL

  49. Manakov SA, Grant SG, Enright AJ: Reciprocal regulation of microRNA and mRNA profiles in neuronal development and synapse formation.

    BMC Genomics 2009, 10:419. OpenURL

  50. Nunez-Iglesias J, Liu CC, Morgan TE, Finch CE, Zhou XJ: Joint genome-wide profiling of miRNA and mRNA expression in Alzheimer's disease cortex reveals altered miRNA regulation.

    PLoS One 5(2):e8898. OpenURL

  51. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS: Human MicroRNA targets.

    PLoS Biol 2004, 2(11):e363. OpenURL

  52. Cheng C, Fu X, Alves P, Gerstein M: mRNA expression profiles show differential regulatory effects of microRNAs between estrogen receptor-positive and estrogen receptor-negative breast cancer.

    Genome Biol 2009, 10(9):R90. OpenURL

  53. van Dongen S, Abreu-Goodger C, Enright AJ: Detecting microRNA binding and siRNA off-target effects from expression data.

    Nat Methods 2008, 5(12):1023-1025. OpenURL

  54. Grimson A, Farh KK, Johnston WK, Garrett-Engele P, Lim LP, Bartel DP: MicroRNA targeting specificity in mammals: determinants beyond seed pairing.

    Mol Cell 2007, 27(1):91-105. OpenURL

  55. Goentoro L, Shoval O, Kirschner MW, Alon U: The incoherent feedforward loop can provide fold-change detection in gene regulation.

    Mol Cell 2009, 36(5):894-899. OpenURL

  56. Tsang J, Zhu J, van Oudenaarden A: MicroRNA-mediated feedback and feedforward loops are recurrent network motifs in mammals.

    Mol Cell 2007, 26(5):753-767. OpenURL

  57. Martinez NJ, Ow MC, Barrasa MI, Hammell M, Sequerra R, Doucette-Stamm L, Roth FP, Ambros VR, Walhout AJ: A C. elegans genome-scale microRNA network contains composite feedback motifs with high flux capacity.

    Genes Dev 2008, 22(18):2535-2549. OpenURL

  58. Tu K, Yu H, Hua YJ, Li YY, Liu L, Xie L, Li YX: Combinatorial network of primary and secondary microRNA-driven regulatory mechanisms.

    Nucleic Acids Res 2009, 37(18):5969-5980. OpenURL

  59. Dai Z, Chen Z, Ye H, Zhou L, Cao L, Wang Y, Peng S, Chen L: Characterization of microRNAs in cephalochordates reveals a correlation between microRNA repertoire homology and morphological similarity in chordate evolution.

    Evol Dev 2009, 11(1):41-49. OpenURL

  60. Chen K, Rajewsky N: The evolution of gene regulation by transcription factors and microRNAs.

    Nat Rev Genet 2007, 8(2):93-103. OpenURL

  61. Herranz H, Cohen SM: MicroRNAs and gene regulatory networks: managing the impact of noise in biological systems.

    Genes Dev 24(13):1339-1344. OpenURL

  62. Hilgers V, Bushati N, Cohen SM: Drosophila microRNAs 263a/b confer robustness during development by protecting nascent sense organs from apoptosis.

    PLoS Biol 2010, 8(6):e1000396. OpenURL

  63. Vergauwen L, Benoot D, Blust R, Knapen D: Long-term warm or cold acclimation elicits a specific transcriptional response and affects energy metabolism in zebrafish.

    Comp Biochem Physiol A Mol Integr Physiol 2010, 157(2):149-157. OpenURL

  64. Logan CA, Somero GN: Transcriptional responses to thermal acclimation in the eurythermal fish Gillichthys mirabilis (Cooper 1864).

    Am J Physiol Regul Integr Comp Physiol 2010, 299(3):R843-852. OpenURL

  65. Kultz D: Molecular and evolutionary basis of the cellular stress response.

    Annu Rev Physiol 2005, 67:225-257. OpenURL

  66. Gualerzi CO, Giuliodori AM, Pon CL: Transcriptional and post-transcriptional control of cold-shock genes.

    J Mol Biol 2003, 331(3):527-539. OpenURL

  67. Yamaguchi-Shinozaki K, Shinozaki K: Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses.

    Annu Rev Plant Biol 2006, 57:781-803. OpenURL

  68. Dresios J, Aschrafi A, Owens GC, Vanderklish PW, Edelman GM, Mauro VP: Cold stress-induced protein Rbm3 binds 60S ribosomal subunits, alters microRNA levels, and enhances global protein synthesis.

    Proc Natl Acad Sci USA 2005, 102(6):1865-1870. OpenURL

  69. Lopez-Maury L, Marguerat S, Bahler J: Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation.

    Nat Rev Genet 2008, 9(8):583-593. OpenURL

  70. Zheng D, Kille P, Feeney GP, Cunningham P, Handy RD, Hogstrand C: Dynamic transcriptomic profiles of zebrafish gills in response to zinc depletion.

    BMC Genomics 2010, 11:548. OpenURL

  71. Zuker M, Stiegler P: Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information.

    Nucleic Acids Res 1981, 9(1):133-148. OpenURL

  72. Shalgi R, Lieber D, Oren M, Pilpel Y: Global and local architecture of the mammalian microRNA-transcription factor regulatory network.

    PLoS Comput Biol 2007, 3(7):e131. OpenURL

  73. Krek A, Grun D, Poy MN, Wolf R, Rosenberg L, Epstein EJ, MacMenamin P, da Piedade I, Gunsalus KC, Stoffel M, et al.: Combinatorial microRNA target predictions.

    Nat Genet 2005, 37(5):495-500. OpenURL

  74. Hammell M, Long D, Zhang L, Lee A, Carmack CS, Han M, Ding Y, Ambros V: mirWIP: microRNA target prediction based on microRNA-containing ribonucleoprotein-enriched transcripts.

    Nat Methods 2008, 5(9):813-819. OpenURL

  75. Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB: Prediction of mammalian microRNA targets.

    Cell 2003, 115(7):787-798. OpenURL

  76. Maslov S, Sneppen K: Specificity and stability in topology of protein networks.

    Science 2002, 296(5569):910-913. OpenURL