Skip to content

Advertisement

You're viewing the new version of our site. Please leave us feedback.

Learn more

BMC Genomics

Open Access

MicroRNAs regulate gene plasticity during cold shock in zebrafish larvae

Contributed equally
BMC Genomics201617:922

https://doi.org/10.1186/s12864-016-3239-4

Received: 17 May 2016

Accepted: 1 November 2016

Published: 15 November 2016

Abstract

Background

MicroRNAs (miRNAs) are critical regulators responding to acute environmental stresses in both plants and animals. By modulating gene expression, miRNAs either restore or reconstitute a new expression program to enhance cell tolerance to stresses. Cold shock is one of the stresses that can induce acute physiological responses and transcriptional changes in aquatic creatures. Previous genomic studies have revealed many cold-affected genes in fish larvae and adults, however, the role of miRNAs in acute cold response is still ambiguous. To elucidate the regulatory roles of miRNAs in the cold-inducible responses, we performed small RNA-seq and RNA-seq analyses and found potential cold regulatory miRNAs and genes. We further investigated their interactions and involvements in cold tolerance.

Results

Small RNA-seq and RNA-seq identified 29 up-/26 down-regulated miRNAs and 908 up-/468 down-regulated genes, respectively, in responding to cold shock for 4 h at 18 °C. miRNA and transcriptomic analyses showed these miRNAs and mRNAs are involved in similar biological processes and pathways. Gene ontology enrichment analyses revealed the cold-induced genes were enriched in pathways, including melanogenesis, GnRH pathway, circadian rhythm, etc. We were particularly interested in the changes in circadian clock genes that affect daily metabolism. The enrichment of circadian clock genes was also observed in previous fish cold acclimation studies, but have not been characterized. To characterize the functional roles of circadian clock genes in cold tolerance, we individually overexpressed selected clock genes in zebrafish larvae and found one of the core clock genes per2 resulted in better recovery from cold shock. In addition, we validated the interaction of per2 with its associate miRNA, dre-mir-29b, which is also cold-inducible. It suggests the transcription of per2 can be modulated by miRNA upon cold shock.

Conclusions

Collectively, our observations suggest that miRNAs are fine turners for regulating genomic plasticity against cold shock. We further showed that the fine tuning of core clock gene per2 via its associated miRNA, dre-mir-29b, can enhance the cold tolerance of zebrafish larvae.

Keywords

Deep sequencingmiRNAClock genesZebrafishEmbryonic developmentCold stress

Background

The survival of aquatic fish is constantly challenged by the changes in environmental factors, such as salinity, oxygen concentration, pH, osmotic pressure and temperature. The change in water temperature has a dramatic effect on physiology and behaviors of aquatic animals that it has been shown to be an “abiotic master factor” for fish [1, 2]. Cold shock is characterized as an acute reduction in ambient temperature that may cause rapid falling in body temperature and a cascade of physiological and behavioral responses [2]. The cold shock stress response in fish can be broadly grouped into three categories; primary, secondary, and tertiary responses [24]. The primary responses of cold shock rapidly change the neuroendocrine pathways, such as catecholamines and corticosteroids [57]. It is initiated by the central nervous system sensing low temperatures that stimulates the chromaffin cells to release catecholamines. The hypothalamus-pituitary-interrenal axis is then activated to secret corticotrophin-releasing hormone to stimulate adrenocorticotrophic hormone and subsequent release of cortisol [4, 8, 9]. Secondary response can be elicited by the primary response. The stress hormones and cortisol interact to induce glucose production in fish via glycogenolysis and glucogenesis [9, 10] that results in changes of protein expression, molecular responses, osmoregulatory, and immune systems [2, 11]. The tertiary responses are the global effects of cold shock to impair development rate, respiratory function, swimming behavior, and mortality [2].

Expressions of numerous genes are regulated by stresses in fish. Transcriptome analysis has been used in the discovery of many stress-inducible genes in cold - acclimated [1214] or cold-shocked fishes [15, 16]. Those genes belong to a variety of biological processes, including chaperones, transcription factors, signal transduction pathways, metabolism, responses to stress or ion transporters like kir2.2 channel in rainbow trout [17], heat shock protein 70 in silver sea bream [18] and glycerol-3-phosphate dehydrogenase in rainbow smelt [19]. However, studies on non-model fishes have been difficult due to the lack of whole genome sequences and inability for genome manipulation. Zebrafish is a favorable model fish due to the availability of whole gene sequences and accessibility of genome manipulation. They can also be used for genomic research using high-throughput RNA-seq and small RNA-seq techniques. Moreover, zebrafish is a tropical animal acclimated to extreme temperature changes ranging from 6.2 to 39.2 °C in certain seasons [20]. This makes it an ideal model to investigate the regulatory mechanisms of temperature-dependent responses.

MicroRNAs (miRNAs) are short non-coding RNAs of around 22 nucleotides, which can modulate the stability and silence mRNAs [21]. A miRNA typically targets multiple mRNAs [22], which makes it a potential fine-tuner for different situations. Studies indicate that miRNAs play a key role in gene regulation in stress responses like oxidative reaction [23], temperature [24] or DNA damage stresses [25] in different organisms [26, 27]. However, relatively less is known about the functions of miRNA in aquatic species. miR-8 and miR-429 have been reported to be involved in the osmotic stress response in zebrafish [28, 29]. In contrast, high-throughput miRNAtomes analysis revealed no significant impact by miRNA-mediated gene regulation in the cold - acclimated adult zebrafish brains [30]. This result is puzzling, but we reasoned that the cold-induced responses might be different in organisms or tissues used. For example, fish larvae are more voluntary to environmental stresses like acute cold shock, but their responses to the changes in miRNA and transcriptomic profiles during acute cold shock remains unexplored.

Yang et al. claimed that miRNAs may have a minor role in regulating transcriptome plasticity in cold acclimated zebrafish brains based on the poor correlations between cold-induced transcriptomic and miRNAtome profiles [30]. This is in clear contrast to our previous understanding observed in other organisms [31, 32]. The discrepancy may be due to differences in organisms, tissue analyzed and experimental conditions. To further clarify this issue, we aimed to investigate the correlation between miRNAtome and transcriptome under cold shock in zebrafish larvae and to resolve possible cold-dependent regulatory mechanisms of selected miRNAs and their target genes. We discovered that similar biological processes and pathways were affected by cold in both miRNAtome and transcriptome analyses. In particular, we found notable correlations between circadian rhythm-related genes (known as circadian clock genes) [33] and their associate miRNAs. Core clock gene per2 plays a major role in regulating circadian rhythm and daily metabolism [34]; we found that dre-mir-29b (dre denotes Danio rerio) can target per2 in zebrafish larvae to possibly regulate cold shock response via modulating the expression of per2. More importantly, we demonstrated that the overexpression of per2 in zebrafish larvae can increase their cold tolerance. These results indicate a pivotal role of miRNA in regulating gene plasticity to counter cold stress.

Methods

Fish maintenance

Zebrafish, Danio rerio, wild-type AB strain, were maintained under a 14 h light/10 h dark cycle at 28.5 °C. Embryos were collected by natural spawning and incubated in 0.3× Danieau’s buffer (1× Danieau’s buffer: 58 mM NaCl, 0.7 mM KCl, 0.4 mM MgSO4, 0.6 mM Ca(NO3)2, and 5.0 mM HEPES in double - distilled water, pH adjusted at 7.6) until examination or fixation.

Cold shock and embryo collection

For cold shock experiment, 50 larvae were cultured at 28.5 °C until 96 h post fertilization (hpf), transferred into a 250 ml beaker containing 100 ml pre-chilled 0.3× Danieau’s buffer, and immediately incubated at designated temperatures until examination. The 96-hpf embryos before cold shock were collected as the control group. Embryos incubated at 28.5 °C or 18 °C for 4 more hours were collected as the normal or cold shock group, respectively. Two different batches of samples were collected for further next generation sequencing (NGS) analysis.

miRNA library construction and sequencing

RNA samples were collected from the control, normal and cold shock groups for cDNA library construction. RNAs of desired time points were isolated by using the TRIzol® Reagent (Life technologies, Carlsbad, CA), and then the amount and quality of RNAs were assessed using the NanoDrop® ND-1000 Spectrophotometer (ThermoFisher, Waltham, MA), Bioanalyzer 2100, RNA 6000 Nano Kit and Small RNA Chip Kit (Agilent Technology, Santa Clara, CA). For miRNA library construction and sequencing, miRNAs were enriched using the PureLink™ miRNA Isolation Kit (Invitrogen, Carlsbad, CA), and then adaptors were hybridized and ligated to the 5′ and 3′ ends of these RNAs using the SOLiD™ Total RNA-Seq Kit (Applied Biosystems Inc., ABI, Carlsbad, CA). The ligated RNAs were used as templates for cDNA synthesis followed by cDNAs purification using the MinElute® PCR Purification Kit (QIAGENE, Valencia, CA). The cDNAs at 120–130 base pair (bp) were selected (Invitrogen, Novex® pre-cast gel products), amplified (ABI, SOLiD™ Total RNA-Seq Kit), purified (Invitrogen, PureLink™ PCR Micro Kit), and proceeded with SOLiD™ System template bead preparation (ABI, SOLiD™ RNA Barcoding Kit). The libraries were sequenced with a read length of 35 bp using the SOLiD 5500 xl sequencer.

cDNA library construction and sequencing

RNA samples were collected as described for cDNA library construction. RNAs of desired time points were isolated as described and used RiboMinus™ (Invitrogen) to remove ribosomal RNAs. For the preparation of cDNA libraries, the RNAs were fragmented using the RNAase III by using the Library Builder™ Whole Transcriptome Core Kit (ABI) and cleaned up by using the Fragmented RNA Concentrator Module Kit (ABI). To amplify and purify the cDNA libraries we used the Library Builder™ Whole Transcriptome Core Kit (ABI). The cDNA libraries with a size between 200 and 300 bp were proceeded with the SOLiD™ System template bead preparation using the EZ Bead™ System (ABI). The libraries were sequenced with read length of 75 bp using the SOLiD 5500 xl sequencer.

Identification of differential expressed genes and miRNAs

Raw reads from the SOLiD 5500 xl sequencing were analyzed using the LifeScope™ Genomic Analysis Software. Those reads were filtered for high quality, adaptor self-ligation, transfer RNAs and ribosomal RNAs (http://rfam.sanger.ac.uk/). These filtering steps resulted in high quality filtered reads representing zebrafish RNA sequences. The filtered reads were mapped to the zebrafish genome sequences (Zv9, http://www.sanger.ac.uk/resources/downloads/zebrafish/) and annotated either the miRBase v18 (http://www.mirbase.org/) for miRNAs or the zebrafish genome annotation file (Zv9/danRer7 from UCSC) for transcripts. Annotated reads were further normalized to the reads per kilo base per million (RPKM). The differentially expressed miRNA and genes were identified and compared using the Partek Genomics Suite (http://www.partek.com/).

Quantitative PCR analysis

The total RNAs were isolated as described and genomic DNAs were removed by using the DNA-Free kit (Life technologies). One microgram of total RNAs was reversely transcribed to miRNAs and cDNAs by the miScript II RT Kit (QIAGENE). The quantitative polymerase chain reaction (qPCR) was performed using the miScript SYBR Green PCR Kit (QIAGENE) for miRNA quantification. The zebrafish U6 (5′-ACTAAAATTGGAACGATACAGAGA-3′) and ef1α served as the internal control for miRNA and mRNA, respectively. qPCR analyses were performed by using the iQ™ SYBR® Green Supermix in a CFX96 station (Bio-Rad, Hercules, California). Primers used are listed in the Additional file 1.

Gene cloning

Coding sequences of period 2 (per2, ENSDARG00000034503), aryl hydrocarbon receptor nuclear translocator-like 1a (arntl1a, ENSDARG00000006791), basic helix-loop-helix family, member e41 (bhlhe41, ENSDARG00000041691) were amplified from the zebrafish total cDNAs by reverse transcriptase PCR (RT-PCR). Total RNAs of zebrafish embryos were isolated as described and reversely transcribed by the MMLV reverse transcriptase (Promega, Fitchburg, Wisconsin). Primers were designed according to the respective reference RNA sequence. The amplicons were sub-cloned into a pT2MUAS vector for Tol2-mediated transgenesis [35]. The vector also contains a viral 2A peptide sequence for the separation of the expressed protein from a nuclear reporter H2AmCherry used for founder fish screening.

Ectopic expression and analysis of circadian clock genes in zebrafish larvae

One-cell stage zebrafish embryos were injected with 25 pg Gal4 mRNAs, 25 pg Toll2 mRNA and 50 pg plasmids and incubated under a 14 h/10 h dark/light cycle until examination. Embryos with H2AmCherry expression were selected for further analysis.

Swimming behavior analysis

Five days post fertilization (dpf) larvae expressing with or without per2, arntl1a or bhlhe41 were individually distributed into each well of a 24-well plate with 1.5 ml of 0.3× Danieau’s buffer and incubated for 4 h at 28.5 or 18 °C. Following incubation, the larvae were transferred to 28.5 °C and videotaped in a controlled light-on DanioVision observation chamber for 30 min. The tracking of swimming for each larva was analyzed during 0–10 min, 10–20 mins and 20–30 mins using the EthoVision XT 8.5 software (Noldus Information Technology, Wageningen, The Netherlands).

Measurement of glucose concentration

Control or transient expressing embryos were cultured and cold shocked as described. Ten larvae for each treatment were collected at 0, 2, and 4 h after cold-shock, lysed in 200 μl RIPA buffer (ThermoFisher Scientific, Waltham, MA), and the glucose concentration was measured by using the AmpliteTM Fluorimetric Glucose Quantitation Kit (AAT-bioquest, Sunnyvale, CA).

Reporter assay of dre-mir-29b

The 3′ untranslated region (UTR) of per2 containing dre-mir-29b target site was amplified from the cDNA library by the forward (F) and reverse (R) primers of the following sequences. F: Bgll2-Stop-per2 ATCGAGATCTTAAAATTCCTTTCGCATTCACAA and R: SacI-per2-ATCGGAGCTCAGTCTGTGAGATCAGTTAAACCA. The amplicons were cloned into the pEGFP-C1 vector by selected restriction enzyme site SacI and Bgll2. Plasmids were mixed with or without dre-mir-29b morpholino (5′ACACTGATTTCAAATGGTGCTAGAT3′). The plasmid mix was injected into 1-cell zebrafish embryos and observed under epifluorescent microscopy at 10 hpf.

Statistical analysis

All experimental values are presented as mean ± standard error and were analyzed by unpaired-sample Student’s t-test and one-way ANOVA.

Results

Cold-induced responses in zebrafish larvae

To investigate cold-induced responses in zebrafish, 96 hpf larvae were exposed to decreasing temperatures to monitor their mobility and mortality for 10 h (Table 1). Fish were mostly normal from 28.5 to 20 °C except for a few inactive fish at 20 °C. In contrast, most fish became motionless after longer incubation at 18 °C. Fish were inactive at 16 °C and were unresponsive upon touching. Death appeared at 12 °C and below. The mortality increased with lower temperatures and longer incubation (data not shown). Thus, we used 18 °C to study cold shock responses to avoid mortality and undesired health problems.
Table 1

Cold impairs swimming and causes mortality in zebrafish larvae

Incubation temperature

Hours Post Treatment

°C

2

4

6

8

10

28.5

Normal

Normal

Normal

Normal

Normal

24

Normal

Normal

Normal

Normal

Normal

20

Normal

Normal

mostly normal but few inactive larvae

Mostly normal but few inactive larvae

Mostly normal but few inactive larvae

18 (% of inactive fish)

83

95

96

100

100

16

All fish were inactive and showed unbalance moves when disturbed

12 (% of motility)

0

0

0

7.2

12

Four - day - old zebrafish larvae were incubated at decreasing temperatures from 28.5 to 12 °C. The mobility and viability were recorded every 2 h until 10 h post treatment

To confirm that the exposure to 18 °C is able to elicit cold responses at the transcriptional level, we examined the change in gene expression of cold - inducible RNA binding protein (cirbp), a well-known cold inducible gene [11]. qPCR analysis showed that the cirbp expression was initially decreased at 2 h post induction (hpi), but significantly bounced back to about 2 folds at 4 hpi and then declined to normal at 10 hpi (Fig. 1a). It appeared that a 4-h incubation at 18 °C is sufficient to induce cold-dependent transcriptomic changes in zebrafish larvae.
Fig. 1

Cold shock treatment in zebrafish larvae. a Relative expression levels of cirbp gene in zebrafish larvae at 96 h post fertilization (hpf) were determined by q PCR at different time point post transfer from 28.5 to 18 °C. * and ** indicates p ≤ 0.05 and 0.01, respectively, comparing to the control at 0 h post transfer. b Scheme of cold exposure. Embryos were incubated at 28.5 °C until 96 hpf then transferred to 18 or 28.5 °C for an additional 4 h. Larvae before (control) and after treatments (normal and cold shock) were then collected for RNA-seq and miRNA-seq analyses as described in the Methods

qPCR analyses detected no notable change in expression of cirbp in samples before and immediately after transferring to 18 °C (data not shown). So the only time 0 control used was samples before cold shock. Larvae were then incubated at 28.5 °C (normal) or 18 °C (cold shock) and samples were taken from both conditions at 4 hpi (Fig. 1b). These experiments were repeated twice and a total of 6 different samples were subjected to small RNA and RNA NGS.

Cold-induced miRNA sequencing profiles

To clarify the roles of miRNA on transcriptomic plasticity in cold responses, we performed both small RNA-seq and RNA-seq analyses in zebrafish larvae. Six RNA samples were extracted from three groups (control, normal and cold shock) as previously described, enriched for 18–30 bp small RNAs for the cDNA library construction and subjected to NGS using the SOLiD 5500 xl. A total of 39,560,685; 17,574,255, and 27,708,966 sequence reads from the first batch were obtained from the control, normal and cold shock sample, respectively (Fig. 2a, the total reads of the second batch were listed in Additional file 2: Figure S2A). The sequencing data was analyzed using the LifeScope™ Genomic Analysis Software and normalized using the Partek Genomics Suite. The length distribution of each library reads were ranged between 18 and 30 bp, and the major sizes peaked around 21~24 bp that matched the general size distribution of miRNAs (Additional file 2: Figure S1A, S2B). Abundant reads were located at chromosome 4, 10, 20, and 25 in different libraries (Additional file 2: Figure S1B, S2C). After filtering out short and adaptor-contaminated reads, about 17–19% miRNAs were identified. The reads and the percentage of different RNA components (miRNA, noncoding RNA, mRNA, novel RNA) were similar among samples (Fig. 2a). To examine the overall effects of cold shock on miRNAtome, the batches of data were normalized by reads count and subjected to principle component analysis to perform the quality check of samples as shown in Additional file 2: Figure S3. Then we mapped the sequencing reads to the miRBase v18, which contains 344 zebrafish miRNAs and found 261 known miRNAs. After normalizing two batches of data from three groups, the reads from different treatments were highly correlated (Fig. 2b). It suggests that the majority of miRNAs remained relatively constant during development and cold shock. Figure 2c shows the log plots of the fold of change in miRNA expression between treatments. We considered the fold of change ≥1.5 folds to be significant. Transcriptomic profiles are changing rapidly during embryonic development. The altered miRNAs between control and normal were thought to be development-dependent. After taking out development-dependent miRNAs, we identified 29 cold-induced miRNAs and 26 cold-repressed miRNAs (within the red boundary in Fig. 2d). The heat map of miRNAs with a change ≥1.5 folds are shown in Fig. 2e and Additional file 3. The top 5 mRNAs affected are listed in Table 2 and the list of affected - miRNAs can be found in Additional file 4.
Fig. 2

Analysis of cold shock - affected miRNAs by small RNA-seq in zebrafish larvae. a A pipeline is presented for small RNA identification at each filtration step (left). The number of reads and percentage passed each filtration step of each treatment are shown on the right chart. b Scatter plots show comparisons of RPKM between treatments. Log plots are shown in (c) for those cold-affected genes with a change in expression ≥2 folds. The vertical line in each graph indicates where log2 (RPKM) equal to 1.07. Only dots with values greater than 1.07 were subjected to further analyses. d The numbers of cold shock-affected (≥1.5 folds changes) miRNA are presented in Venn diagram as shown in Fig. 1e. e Heat maps profiles for change in expression of selected miRNAs

Table 2

Cold shock-induced differential expression of known miRNAs

miRNA

Mature miRNA sequence

Fold change

Up-regulated miRNAs

 Dre-mir-29b-1

UAGCACCAUUUGAAAUCAGUGU

3.90

 Dre-mir-99-2

AACCCGUAGAUCCGAUCUUGUG

2.04

 Dre-mir-99-1

AACCCGUAGAUCCGAUCUUGUG

1.79

 Dre-mir-92a-2

UAUUGCACUUGUCCCGGCCUGU

1.75

 Dre-mir-2184

AACAGUAAGAGUUUAUGUGCU

1.73

Down-regulated miRNAs

 Dre-mir-737

AAUCAAAACCUAAAGAAAAUA

−1.71

 Dre-mir-9-2

UCUUUGGUUAUCUAGCUGUAUGA

−1.75

 Dre-mir-363

AAUUGCACGGUAUCCAUCUGUA

−1.79

 Dre-mir-125b-2

UCCCUGAGACCCUAACUUGUGA

−1.91

 Dre-mir-199-1

CCCAGUGUUCAGACUACCUGUUC

−2.79

After eliminating miRNAs whose expression might be changed due to development process during incubation (control vs normal). The top 5 of 29 up- and 26 down-regulated known miRNAs (annotated in miRBase v18) are shown, along with their sequences and fold changes compared to the control groups

Validation of miRNA-seq data with qPCR

We performed qPCR analysis to validate the changes of the cold-affected miRNAs (dre-mir-29b, the most up-regulated miRNA, and 9 other randomly selected miRNAs) and compared them to those observed in small RNA-seq results (Fig. 3a). The Pearson’s correlation coefficient between two analyses is 0.87. This is relatively high and suggests the reliability of small RNA-seq results (Fig. 3b).
Fig. 3

Validation of small RNA-seq data and biological process pathway summary for the target genes of significantly affected miRNA upon cold shock. a Comparison of change in relative expression of selected miRNAs determined by NGS and qPCR (N = 3) (b) Pearson’s correlation between small RNA-seq and qPCR data (correlation coefficient (R2 = 0.807). c Pie charts show the proportions of cold-induced (left) and repressed (right) miRNAs whose target genes are involved in different biological processes with a p value ≤ 0.05

Function annotation for target genes of the cold-dependent miRNAs

The cold-dependent target genes of each miRNA were predicted by TargetScanFish release 6.2 (http://www.targetscan.org/fish_62/) with the total context+ score −0.3 as a cutoff [36]. We obtained 6263 and 5381 targets of cold-induced and -repressed miRNAs, respectively (Additional file 4). We then conducted the target genes functional annotation using The Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/) [37]. The top 10 biological processes were ranked according to P value (Table 3 and Additional file 5). The P value is an indication of gene enrichment in annotated biological process. The percentages of affected miRNA target genes with a value ≤ 0.005 are illustrated in pie graphs (Fig. 3c). In general, both up- and down-regulated miRNA target genes were enriched in processes like metabolism, transport, phosphorylation, development, signaling transduction, sexual determination and reproduction. In contrast, transcription-related genes were only enriched in the up-regulated miRNA target genes, but cell movement and motility-related genes were abundant in the down-regulated miRNA target genes.
Table 3

The top 10 biological processes for cold-related miRNA target genes

Biological Processes of

P value*

Cold-induced miRNA target genes

 Regulation of transcription, DNA-dependent

3.30E-04

 Polyol metabolic process

3.40E-04

 Regulation of RNA metabolic process

5.20E-04

 Regulation of transcription

6.70E-04

 Regulation of small GTPase mediated signal transduction

1.10E-03

 Ion transport

1.10E-03

 Glycerol metabolic process

1.30E-03

 Alditol metabolic process

1.30E-03

 Organophosphate metabolic process

1.50E-03

 Metal ion transport

1.60E-03

Cold-repressed miRNA target genes

 Metal ion transport

1.60E-04

 Protein amino acid phosphorylation

3.60E-04

 Intracellular signaling cascade

7.60E-04

 Ion transport

8.50E-04

 Cation transport

1.80E-03

 Potassium ion transport

1.80E-03

 Regulation of small GTPase mediated signal transduction

1.80E-03

 Organophosphate metabolic process

3.10E-03

 Monovalent inorganic cation transport

3.60E-03

 Regulation of Ras protein signal transduction

4.70E-03

*calculated by Fisher Exact which is adopted to measure the gene-enrichment in annotation biological process

Target genes of 29 up- and 26 down-regulated known miRNAs were predicted by TargetScan (Additional file 5). The target genes were annotated by DAVID bioinformatics resources 6.7. The top 10 biological processes of the target genes and P values are shown

We also predicted the top 10 cold-induced and cold-repressed miRNA-mediated signaling pathways by using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Table 4 and Additional file 5). In the up-regulated group, most genes were involved in metabolism, including cysteine and methionine metabolism, fructose and mannose metabolism, glycerophospholipid metabolism, inositol phosphate metabolism and sphingolipid metabolism. One pathway was related to extracellular matrix-receptor interaction and the other one was related to biosynthesis of unsaturated fatty acids. In the cold-repressed group, many target genes were involved in biosynthesis and metabolism pathways, including chondroitin sulfate biosynthesis, synthesis and degradation of ketone bodies, O-Glycan biosynthesis, glycosphingolipid biosynthesis, biosynthesis of unsaturated fatty acids, sphingolipid metabolism and alanine, aspartate and glutamate metabolism. Some cold-repressed miRNA target genes were also found in the phosphatidylinositol, notch and lysozyme signaling pathways.
Table 4

The top 10 KEGG pathways for cold-related miRNA target genes

KEGG pathway of

P value*

Cold-induced miRNA target genes

 Cysteine and methionine metabolism

0.0027

 ECM-receptor interaction

0.0035

 Fructose and mannose metabolism

0.005

 Biosynthesis of unsaturated fatty acids

0.011

 Glycerophospholipid metabolism

0.011

 Inositol phosphate metabolism

0.014

 Sphingolipid metabolism

0.022

 Melanogenesis

0.028

 GnRH signaling pathway

0.029

 Circadian rhythm

0.034

Cold-repressed miRNA target genes

 Chondroitin sulfate biosynthesis

0.00081

 Synthesis and degradation of ketone bodies

0.0092

 Phosphatidylinositol signaling system

0.017

 Sphingolipid metabolism

0.019

 Notch signaling pathway

0.023

 O-Glycan biosynthesis

0.025

 Lysosome

0.03

 Alanine, aspartate and glutamate metabolism

0.035

 Glycosphingolipid biosynthesis

0.036

 Biosynthesis of unsaturated fatty acids

0.038

*calculated by Fisher Exact which is adopted to measure the gene-enrichment in annotation KEGG pathway

Target genes of 29 up- and 26 down-regulated known miRNAs were predicted by TargetScan (Additional file 5). The targets genes were annotated by DAVID bioinformatics resources 6.7. The top 10 KEGG pathway of the target genes and P values are shown

Cold-induced RNA sequencing profile

To validate our predicted miRNA targets, six RNA samples were extracted from zebrafish larvae, treated and processed as described previously to perform RNA-seq analysis using the SOLiD 5500 xl. The transcriptomic differences between three groups of fish were similarly analyzed as described in the miRNAtome analysis. Total numbers of reads and mapped reads of duplicates from three different treatments are shown in Fig. 4a. After normalizing two batches of data from three groups, the RPKM from different treatments (Fig. 4b) were highly correlated. It suggests that the majority of mRNAs remained relatively constant during development and cold shock. Figure 4c shows the log plots of the fold of change in RNA expression between treatments. We set the fold change ≥ 2 folds and RPKM ≥ 0.1 to be significant. After taking out development-dependent RNAs, we obtained 908 cold-induced genes and 463 cold-repressed genes (within the red boundary in Fig. 4d). The list of affected genes can be found in Additional file 6. The heat map of genes with a change ≥ 2 folds are shown in Fig. 4e and Additional file 3. Among them, the genes related to circadian rhythm are indicated.
Fig. 4

Transcriptomic analysis of cold shock-induced genes in zebrafish larvae. a The number of reads and percentage passed each filtration step of each treatment are shown on the right chart. b Scatter plots show comparisons of reads per kilo base per million (RPKM) between treatments. Log plots are shown in (c) for those cold-affected genes with a change in expression ≥ 2 folds and RPKM ≥ 0.1. The total dot numbers are shown on the right top corner of each graph. d Venn diagrams show the numbers of cold shock-induced and repressed genes between treatments. e Heat maps profiles for mRNA expression of selected genes with annotations of circadian - related genes

Validation of the cold-induced changes in circadian clock genes with qPCR

The circadian clock genes were also shown to be affected during cold acclimation in zebrafish [38]. However, their roles in the cold-induced response have not been elucidated. We first performed qPCR analysis to validate the changes of those circadian clock genes, which include per2, per1a, clock3, clock, arntl1a, arntl1b, nr1d, cry1a, and bhlhe41, under cold shock or normal conditions. The qPCR data was compared to those results by the RNA-seq (Fig. 4a). The Pearson’s correlation coefficient between two analyses is 0.807 (Fig. 4b) that is relatively high and suggests the reliability of RNA-seq results.

Function annotation of cold-dependent genes in transcriptome

To reveal the biological interpretation of cold-dependent genes, we conducted the genes’ functional annotation using DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/) to find the enriched biological processes by the Gene Ontology (GO) term analysis. The P value ranking indicates gene enrichment in annotated biological processes. The top 10 biological processes, such as regulation of RNA metabolic process, are listed in Table 5 and the number of genes belonging to each process are shown in Fig. 5c. The percentages of affected genes with a value ≤ 0.1 are illustrated in a pie graph (Fig. 5d). In general, both up- and down-regulated genes were enriched in regulation of RNA metabolic processes, regulation of transcription, DNA-dependent and regulation of transcription. In contrast, protein amino acid phosphorylation, phosphate metabolic process, calcium ion transport, neuron differentiation and pattern specification process were only enriched in the up-regulated genes. Response to DNA damage stimulus, cellular response to stress, cell fate commitment, embryonic morphogenesis, inorganic anion transport, anion transport and DNA repair were abundant in the down-regulated genes.
Table 5

The top 10 biological processes for cold-related genes

Biological Processes of

P value*

Cold-induced genes

 Regulation of RNA metabolic process

1.50E-03

 Regulation of transcription, DNA-dependent

1.80E-03

 Protein amino acid phosphorylation

1.90E-03

 Phosphate metabolic process

2.10E-03

 Rgulation of transcription

5.40E-03

 Phosphorylation

6.80E-03

 Transcription

9.50E-03

 Calcium ion transport

0.01

 Neuron differentiation

0.015

 Pattern specification process

0.019

Cold-repressed genes

 Response to DNA damage stimulus

0.01

 Cellular response to stress

0.018

 Regulation of transcription

0.019

 Cell fate commitment

0.027

 Regulation of RNA metabolic process

0.03

 Embryonic morphogenesis

0.033

 Inorganic anion transport

0.037

 Regulation of transcription, DNA-dependent

0.04

 Anion transport

0.04

 DNA repair

0.047

*calculated by Fisher Exact which is adopted to measure the gene-enrichment in annotation biological process

The cold-related genes (Additional file 6) were annotated by DAVID bioinformatics resources 6.7. The top 10 biological processes of the cold-related genes and P values are shown

Fig. 5

Validation of RNA-seq data and biological process pathways summary of the genes significantly changed upon cold shock. a Comparison of change in relative expression of circadian clock genes determined by NGS and qPCR (N = 3). b Pearson’s correlation result between RNA-seq and qPCR data. (Correlation coefficient = 0.807) (c) Numbers of cold-induced (left) and repressed (right) genes belong to the top 10 biological processes. d Pie charts show the proportions of cold-induced (left) and repressed (right) genes which involve in different biological processes with a p value ≤ 0.05

We also analyzed the cold-affected genes in the KEGG pathway database (Table 6) [39]. In the up-regulated groups, most of the genes were involved in MAPK, hedgehog, GnRH and insulin signaling pathways. Other genes were enriched in immunity-related pathways like retinoic acid-inducible gene 1 (RIG-I)-like receptor and oligomerization domain (NOD)-like receptor, circadian rhythm, melanogenesis and adhesion molecule. In the cold-repressed genes, only neuroactive ligand-receptor interaction and steroid biosynthesis were significantly enriched.
Table 6

The top 10 KEGG pathways for cold-related genes

KEGG pathway of

P value*

Cold-induced genes

 MAPK signaling pathway

2.60E-05

 Circadian rhythm

2.60E-04

 RIG-I-like receptor signaling pathway

0.011

 GnRH signaling pathway

0.025

 Cell adhesion molecules (CAMs)

0.026

 Adherens junction

0.034

 Melanogenesis

0.036

 Hedgehog signaling pathway

0.049

 NOD-like receptor signaling pathway

0.056

 Insulin signaling pathway

0.057

Cold-repressed genes

 Neuroactive ligand-receptor interaction

3.90E-03

 Steroid biosynthesis

6.50E-03

*calculated by Fisher Exact which is adopted to measure the gene-enrichment in annotation KEGG pathway

The cold-related genes (Additional file 6) were annotated by DAVID bioinformatics resources 6.7. The top 10 KEGG pathway of the genes and P values are shown

Correlation of expression profiles of miRNA and mRNA

To determine the correlation between cold-affected miRNAtome and transcriptome, we further analyzed the expression of cold-affected miRNAs and their target genes as shown in Additional file 7. The affected miRNAs should be negatively correlated with their target mRNAs. However, it is not always the case. Only 23.8% of up-regulated miRNAs resulted in down-regulation of target mRNAs. On the other hand, 75.4% of down-regulated miRNAs caused up-regulation of target mRNAs. It appears that the release of miRNA inhibition may have a bigger impact on the expression of target mRNAs.

Overexpression of per2 assists in cold recovery in zebrafish larvae

Biological processes of the Gene Ontology (GO) term analysis shows that circadian rhythm is enriched at the first category - regulation of RNA metabolic process as shown in Additional file 8. As described previously, the roles of circadian clock genes in the cold shock responses are still illusive. In the molecular level, the core clock genes, CLOCK and NPAS2/BMAL1, induce expression of repressors, PER1-3/CRY1-2 to form an interacting feedback loop to sustain the daily oscillation [40]. BHLHE41 plays a relative minor role that functionally resembles the inhibitor of the core clock gene complex and oscillates similarly with the PER1-3/CRY1-2 complex [41]. It has been reported that the dynamics of circadian clock genes are highly correlated with cell metabolism [42, 43]. We hypothesized that the circadian gene-dependent changes in cell metabolism may help in cold resistance. To test this hypothesis, we selected three circadian clock genes, arntl1a (bmal), bhlhe41 and per2, which had different levels of change (1.14, 2.4 and 12.6 folds, respectively) upon cold shock, to perform functional assays. We designed overexpressing plasmids containing Tol2, 5XUSA element, different circadian clock gene coding region sequence and H2AmCherry indicator. The circadian clock gene and mCherry were intervened with a 2A peptide sequence for separation upon translation. Different plasmids were separately injected into 1-cell embryos with gal4 mRNAs (Fig. 6a). The ectopic expression of these plasmids were screened by the H2AmCherry fluorescence under epifluorescent microscope as an indirect indication of circadian clock gene overexpression (Fig. 6b). To validate this, we performed qPCR analyses and found that the expressions of bmal, per2 and bhlhe41 were increased to 18.86 ± 9.1, 12.85 ± 11.17 and 5.02 ± 1.88 folds, respectively (Fig. 6c, N = 3). These fish with transient expression of circadian clock genes were then used for the cold tolerance assay.
Fig. 6

Overexpression of core clock gene per2 or bhlhe41 increases cold tolerance. a The design of the Gal-UAS driven circadian clock gene expression construct is shown. H2AmCherry (red) is placed upstream of clock gene (green) that will be cleaved by a 2A peptide (yellow). b 1-Cell embryos were injected with indicated plasmids and examined under bright (upper row) and dark filed (bottom row through a rhodamine filter cube). Mosaic expression of mcherry indicates the expression of bmal, per2, or bhlhe41 in 48 hpf embryos. c Normalized expression level of clock genes compared to control embryos at 5 days post fertilization by qPCR (N = 3). d Scheme of cold recovery assay at 5 days post fertilization fish. e Classifications of swimming track recorded for 10 mins. f Quantitative analysis of swimming patterns for zebrafish larva transient expressed without (control) and with designated circadian clock gene at different time periods after recording. g Glucose concentration were measured in 4 dpf zebrafish larval lysate collected at different time point after cold shock. *P ≤ 0.05, **p ≤ 0.01

Zebrafish swimming behavior has been used as an indicator for cold tolerance [44]. We first cultured the circadian clock gene overexpressing fish under a 14 h/10 h dark light cycle for 5 days and incubated in 18 °C for 4 h. We then transferred fish to a 28.5 °C observation chamber to record the spontaneous swimming behavior from 0 to 30 mins after transfer (Fig. 6d). The moving paths at different intervals after cold recovery (0–10, 10–20 and 20–30 min) were recorded and analyzed according to the intensities of four classes of swimming activity (Fig. 6e). The per2 and bhlhe41 overexpressing fish had higher swimming activity than those of control and arntl1a overexpressing fish at 10–20 mins interval (% of moving fish: Ctrl: 23 ± 7%, n = 40; Bmal: 35 ± 7%, n = 34; Per2: 39 ± 15%, n = 35, Bhlhe41: 42 ± 15%, n = 36). The per2-overexpressing fish exhibited the most active swimming behavior during the 20–30 min interval (% of moving fish: Ctrl: 44 ± 9.2%, n = 32; Bmal: 42 ± 12%, n = 36; Per2: 76 ± 12%, n = 30; Bhlhe41: 53 ± 10%, n = 34) as shown in Fig. 6f.

To examine cell metabolism, we measured the change in glucose concentration in the circadian clock gene overexpressing larvae. In control and bmal (arntl1a)-overexpressing fish, the glucose concentration was gradually increased during cold shock. In contrast, the bhlhe41- and per2-overexpressing fish had a higher basal glucose level and the level was not significantly changed during cold shock (Fig. 6g).

Dre-mir-29b targets per2 mRNA and regulates its expression dynamically during cold treatment

The increase in cold tolerance in per2-overexpressing fish suggests a pivotal role of circadian clock genes against cold stress. Comparing the expression level of miRNA-seq and mRNA-seq of gene of interest, we identified several cold-dependent miRNAs that may regulate circadian clock genes (Tables 7 and 8). We selected dre-miR-29b to illuminate how miRNA affect transcriptome plasticity during cold shock (Fig. 7a). Dre-miR-29b is the top one cold-induced miRNA with a predicted target of per2. To address this, we injected pEGFPC1-per2 3′UTR with or without dre-mir-29b morpholino oligonucleotides (MO) (Fig. 7b). The co-injection of mir-29b MO significantly enhanced the % of embryos showing the fluorescence from the pEGFPC1-per2 3′UTR (51.9%, n = 95, N = 3) compared to that of plasmid only embryos (8.65%, n = 49, N = 3) as shown in Fig. 7c. This data suggests that dre-mir-29b MO can prohibit dre-mir-29b targeting per2 mRNAs.
Table 7

Cold-affected circadian clock genes and associated miRNAs

Circadian rhythm

miRNA

Basic helix-loop-helix domain containing, class B, 3 like

Dre-mir-92, dre-mir-107, dre-mir-103

Cryptochrome 3 (cry3)

Dre-mir-29

Aryl hydrocarbon receptor nuclear translocator-like 1b (arntl1b)

Dre-mir-99

Clock homolog 3 (clock3)

Dre-mir-29, dre-mir-96, dre-mir-499

Period homolog 2 (per2)

Dre-mir-29,dre-mir-18

Period homolog 1a (per1a)

Dre-mir-135

Nuclear receptor subfamily 1 (nr1d1)

Dre-mir-204

Basic helix-loop-helix family, member e40 (bhlhe40)

Dre-mir-128, dre-mir-181

Casein kinase 1, epsilon

Dre-mir-130, dre-mir-301

Casein kinase 1, delta a

Dre-mir-29

Table 8

Cold-affected miRNAs and corresponding circadian clock genes

MicroRNA

Up/Down

Target

Dre-mir-29b

Up

Per2, Per3, Cry3, Ck1,Clock3

Dre-mir-135

Up

Per1a, Nr1d4a

Dre-mir-99

Up

Arntl1b

Dre-mir-92a

Up

Bhlhe41

Dre-mir-96

Up

Clock3

Dre-mir-204

Up

Nr1d1

Dre-mir-128

Up

Bhlhe40

Dre-mir-107

Up

Bhlhe41

Dre-mir-2192

Down

Clock

Dre-mir-18

Down

Per2, Per3,Bhlhe41

Dre-mir-202

Down

Nr1d2b

Dre-mir-181

Down

Bhlhe40, Nr1d12a

Dre-mir-363

Down

Nr1d2a, Bhlhe41

Dre-mir-130

Down

Ck1ε, Nr1d12a

Dre-mir-301

Down

Ck1ε

Dre-mir-184

Down

Cry1b

Dre-mir-499

Down

Clock3, Nr1d14b

Up/Down: Significantly up-regulated or down-regulated miRNAs are listed, the comprehensive Targeting clock genes are listed at the right column

Fig. 7

Dre-mir-29b may regulate per2 expression during cold acclimation. a per2 is a predicted target of dre-mir-29b. A fragment of 3′ untranslated region (3′UTR) containing the complementary nucleotides of dre-mir-29b is shown. b 1-cell embryos were injected 50 pg pEGFPC1 bearing 3′UTR of per2 target site in the presence or absence of 5 ng mir-29b MO and examined at 90% epiboly stage under epifluorescent microscopy using a FITC cube. c The percentages of EGFP - positive embryos injected without or with dre-mir-29b are shown. The numbers of embryos used in each treatment are presented at the bottom of each column. (N = 3), **p ≤ 0.01. d Embryos were treated with cold shock, lysed at designated time post cold shock, and their RNAs were extracted for qPCR analysis for dre-mir-29b (open bar) and per2 (filled bar)

To further study the dynamic expression of per2 and dre-mir-29b during cold shock, we measured their expression levels at a shorter duration (every 2–4 h) until 24 hpi (Fig. 7d). It appeared that both dre-mir-29b and per2 expression was elevated at 4 hpi. The per2 expression declined afterward. In contrast, the dre-mir-29b expression increased afterwards, but a sudden drop was observed at 24 hpi. These results suggested that dre-mir-29b may act to control the per2 expression-induced by cold shock.

Discussion

Cold shock is a general, but harmful, environmental stress for ectothermic animals. We performed high throughput small RNA-seq and RNA-seq experiments to investigate potential regulatory mechanisms for cold-inducible responses and to characterize cold-affected miRNAs and genes in zebrafish larvae. We found 29 cold-induced miRNAs and 26 cold-repressed miRNAs in the microRNAtome analysis. We also discovered 908 cold-induced genes and 468 cold-repressed genes in the transcriptome analysis. Among them, circadian clock genes and others were enriched in both analyses. Thus, we overexpressed one of the cold-induced core clock gene per2 and showed a notable increase in cold tolerance in zebrafish larvae. We further validated the targeting of per2 by its associated miRNA, dre-mir-29b, which was also induced by cold shock. This suggests that miRNA can serve as a fine tuner for regulating the transcription of cold-inducible genes during cold stress.

An acute temperature drop may result in lethality in aquatic animals. Therefore, we first determined a permissible temperature (18 °C) that could induce adaptive cold responses without deleterious effects on fish. Cirbp, a cold shock-inducible protein, is served as a molecular indicator for cold shock. Cirbp acts as a RNA chaperone to unfold RNA secondary structures to facilitate translation. The elevation of cirbp expression has been known to be a universal marker upon cold exposure [1]. We measured the cirbp expression by qPCR and found that 4 h is the minimum duration which could cause significant gene induction. Furthermore, transcriptome and miRNAtome profiles change rapidly during the larvae stage, however, it was generally ignored in previous studies [15, 38]. Therefore, we ruled out the genes, which might be changed during development to be sure that the up- and down- regulated genes selected are truly cold shock induced.

The cold-affected genes and miRNAs are involved in many crucial biological processes. The most dominate pathways are attributed to transcription regulation, protein phosphorylation and transport. According to previous studies, both cold acclimation and cold shock could affect RNA transcription [11, 15, 45]. In our RNA-seq data, around 20% cold-regulated genes encode transcription factors such as her, hsf2, hsf1, dmrt3a, pou1f1, atf3, gata1a, cebpb, fosl1, bhlhe41, bhlhe40, nr1d4b, nr1d4a, nr1d1 and nr1d2a. Among these genes, bhlhe41, bhlhe40, nr1d4b, nr1d4a, nr1d1, nr1d2a have been reported to be enriched in cold-induced carp larvae [11], suggesting transcription of certain genes should be induced by the cold in the fish larvae stage. The protein phosphorylation landscape is changed during cold treatment in aquatic creatures [46]. In the same token, we also observed the up-regulation of many kinase proteins such as MAP kinase interacting serine/threonine kinase 1 cascades, calcium/calmodulin-dependent protein kinase, protein kinase C and casein kinase 1. Furthermore, we found the changes in expression of multiple ion transporters (atp2a2a, atp2a2b, atp2b3a, cacnb2a, cacna1c, ryr1b, and slc8a2b). It implies that cold shock may inhibit ion pumps and channels that results in imbalance of osmoregulation in fish [47, 48]. Studies also revealed that cold can affect the oxidation reduction processes [15, 45]. However, these processes were not significantly changed in our analyses, which may reflect different cold treatment regime or animals used. In general, our results are in good agreement with previous reports.

Circadian clock genes, including bhlhe40, npas2, nr1d1, nr1d4, per1a, per1b and per2, are under the regulation of RNA metabolic process. The RNA metabolic process is the No. 3 and No. 1 in the cold-induced miRNA target genes (Table 3) and the cold-induced genes (Table 5), respectively. We analyzed genes belonging to the regulation of RNA metabolic process by using the KEGG pathway analysis and found that circadian rhythm is one of two most enriched pathways as shown in the second worksheet of Additional file 8. It gains further proof that circadian rhythm is listed as the No. 2 pathway for all cold-induced genes analyzed as shown in Table 6. Collectively, these results indicate circadian clock genes are notably affected by cold shock.

A group of circadian clock genes, and their associated miRNAs were upregulated upon cold shock (Tables 7 and 8). Similar upregulation of circadian clock genes by cold stress have also been observed in previous studies [15, 38]. It suggests a possible involvement of circadian associated regulation during cold acclimation and cold shock. However, direct evidence for a role of circadian clock genes in cold responses is still lacking. Circadian rhythm controls the daily rhythm for an organism in behavior, physiology and biochemistry [40]. The molecular basis of this oscillation consists of interaction between positive and negative feedback loops. At the molecular level, the circadian rhythm is governed by a set of transcription activators, CLOCK and NPAS2/BMAL1, which induce expression of repressors, PER1-3/CRY1-2, to form an interacting feedback loop [49]. BHLHE41 functions to inhibit CLOCK:BMAL1 transactivation of the clock gene Per1 to resembles the negative feedback loop of PER/CRY complex [41]. In addition to light, circadian rhythm could also be entrained by temperature and feeding behavior [50]. Zebrafish core clock genes, including cry2a, cry3 and per4, can be entrained by lower temperature but not per2 [51]. This implies the elevation of per2 in our study may play a novel biological function in the cold stress response. Here, we successfully demonstrate that overexpression of core clock gene per2 can increase cold tolerance in zebrafish larvae upon cold shock, but not bmal. Previous studies showed that knockout per2 impairs gluconeogenesis [52] but bmal null mice increased gluconeogenesis [53], suggesting two proteins may play opposite role in glucose metabolism. The per2 gene may have a glucocorticoid response element that responds to glucocorticoid hormones to increase blood glucose [54]. The elevation in plasma glucose is enhanced by the glucocorticoid/cortisol-mediated gluconeogenesis, which is a common stress response in aquatic creatures in response to cold [55]. We observed the basal glucose level remains high in transgenic zebrafish larvae overexpressing per2, which suggests that the elevation of per2 expression level may counter cold stress by maintaining a higher gluconeogenesis. Per2 is also known to control lipid metabolism by direct regulation of PPARγ [56]. Whether Per2 also exerts its cold-defending function via enhanced lipid metabolism still awaits further investigation. Furthermore, bhlhe41-overexpressing fish also showed better cold tolerance in locomotor activity and retain higher glucose level. It will be interesting to study whether the inhibition of CLOCK/BMAL complex would enhance cold tolerance.

miRNAs play a key role in cellular responses to the changing environment. Upon stress, cells also tend to restore or reprogram their transcriptomes to counter the changes modulated by miRNAs [27]. It implies a plausible role of miRNAs in mediating cold responses. However, Yang et al. previously reported that cold acclimation trigger changes in some miRNA profiles in adult zebrafish brains and those miRNAs have poor associations to the affected genes. They concluded that microRNAs play a minor role in mediating transcriptome plasticity during cold stress [30]. It appears to be contradictory to our finding that the affected miRNAs have predicted target genes involved in cold-inducible biological processes, including melanogenesis, GnRH signaling pathway and circadian rhythm.

To clarify the interaction between per2 and associated dre-mir-29b, we examined their dynamic expression during cold shock and found good correlations for dre-mir-29b to serve as an inhibitor of per2 expression. The dynamic pattern between miRNAs and target genes nicely matches the “pulse model” as described by Leung et al. [27]. As miRNAs accumulate over time, miRNAs exert their suppressive effect on the mRNA targets, causing a decrease in mRNA target expression. The poor correlation between cold-affected miRNAs and mRNAs observed previously [30] might be due to the time of sample collection or the use of different tissues (brain tissue vs. whole larvae used here). The duration of stress also determines if cells can re-establish homeostasis or initiate the mechanism to acclimate to the new environment [57]. The discrepancy between two studies may also be resulted from different machineries used in two different stress levels (cold acclimation v.s. cold shock).

Conclusions

Here, we demonstrate for the first time that miRNAs can be fine turners for regulating genomic plasticity against cold shock. We further showed that the fine tuning of core clock gene per2 via its associated miRNA, dre-mir-29b, can enhance the cold tolerance of zebrafish larvae. The circadian clock genes may thus serve as targets for designing novel strategies to face cold challenges in aquaculture.

Abbreviations

3′ UTR: 

3′ untranslated region

arntl1a

Aryl hydrocarbon receptor nuclear translocator-like 1a

bhlhe41

Basic helix-loop-helix family, member e41

cirbp

Cold inducible RNA binding protein

DAVID: 

The Database for Annotation, Visualization and Integrated Discovery

Dpf: 

Days post fertilization

Dre: 

Danio rerio

F: 

Forward

GnRH: 

Gonadotropin releasing hormone

GO: 

Gene ontology

H2A: 

Histone 2A

hpf: 

h post fertilization

hpi: 

h post induction

KEGG: 

Kyoto Encyclopedia of Genes and Genomes

miRNAs: 

MicroRNAs

MO: 

Morpholino oligonucleotides

NOD: 

Oligomerization domain

Nts: 

Nucleotides

Per2: 

Period 2

qPCR: 

Quantitative polymerase chain reaction

R: 

Reverse

RIG-I: 

Retinoic acid-inducible gene 1

RPKM: 

Reads per kilo base per million

UTR: 

Untranslated region

Declarations

Acknowledgements

We thank Taiwan Zebrafish core facility for providing technical guidance and plasmids. We thank Technology Commons, College of Life Science, Center for Systems Biology, National Taiwan University for excellent technical assistance with NGS and qPCR experiments. We thank the Taiwan Bioinformatics Core at the National Cheng Kung University, supported by the Taiwan Ministry of Science for assisting with computational analyses and data mining. We thank the Center for Genomic Medicine, National Cheng Kung University, for the support of NGS analyses. The authors would also like to express great appreciation to the staffs of the zebrafish Core at National Taiwan University for providing assistance in fish maintenance. We are also grateful to Dr. Shih-Kuo Chen for the helpful discussion in the regulation of circadian rhythm.

Funding

This work was supported by the Ministry of Science and Technology, Taiwan (101-2321-B-002-054-, 102-2121-B-002-041, 103-2321-B-002-009) and National Taiwan University (NTU NTU CESRP-10R70602A5 and NTU ERP-10R80600 to SJL

Availability of data and materials

No extra dataset except those presented in this article and the supporting information.

Authors’ contributions

ICH carried out the clock gene functional studies, participated in the data analysis and drafted the manuscript. YCH carried out cold shock scheme analysis, miRNAome and transcriptome analysis and drafted the manuscript. HSS designed and supervised the next generation sequence experiments. TMC designed, analyzed and supervised miRNAome and transcriptome experiments, and edited the manuscript. SJL designed, analyzed, supervised all experiments, and revised the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

All experimental procedures on zebrafish were approved by the use of laboratory animal committee at National Taiwan University, Taipei, Taiwan (IACUC Approval ID: 100 Animal Use document No. 92) and carried out in accordance with the approved guidelines.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Department of Life Science, National Taiwan University
(2)
Institute of Molecular Medicine, College of Medicine, National Cheng Kung University
(3)
Institute of Basic Medical Sciences, College of Medicine, National Cheng Kung University
(4)
Department and Graduate Institute of Aquaculture, National Kaohsiung Marine University
(5)
Research Center for Developmental Biology and Regenerative Medicine, National Taiwan University
(6)
Center for Biotechnology, National Taiwan University
(7)
Center for Systems Biology

References

  1. Brett JR. Energetic responses of salmon to temperature. A study of some thermal relations in the physiology and freshwater ecology of sockeye salmon (Oncorhynchus nerka). Am Zool. 1971;11(1):99–113.View ArticleGoogle Scholar
  2. Donaldson MR, Cooke SJ, Patterson DA, Macdonald JS. Cold shock and fish. J Fish Biol. 2008;73(7):1491–530.View ArticleGoogle Scholar
  3. Mazeaud MM, Mazeaud F, Donaldson EM. Primary and secondary effects of stress in fish - some new data with a general review. T Am Fish Soc. 1977;106(3):201–12.View ArticleGoogle Scholar
  4. Barton BA. Stress in fishes: a diversity of responses with particular reference to changes in circulating corticosteroids. Integr Comp Biol. 2002;42(3):517–25.View ArticlePubMedGoogle Scholar
  5. Barton BA. Physiological changes in fish from stress in aquaculture with emphasis on the response and effects of corticosteroids. 1991.Google Scholar
  6. Flos R, Reig L, Torres P, Tort L. Primary and secondary stress responses to grading and hauling in rainbow-trout, salmo-gairdneri. Aquaculture. 1988;71(1-2):99–106.View ArticleGoogle Scholar
  7. Chen WH, Sun LT, Tsai CL, Song YL, Chang CF. Cold-stress induced the modulation of catecholamines, cortisol, immunoglobulin M, and leukocyte phagocytosis in tilapia. Gen Comp Endocrinol. 2002;126(1):90–100.View ArticlePubMedGoogle Scholar
  8. Mazeaud, MM and Mazeaud, F. Adrenergic responses to stress in fish. In: Pickering, AD, editor. Stress and fish. New York: Academic Press; 1981. p. 49-75. Google Scholar
  9. Marcel Martinez-Porchas LRM-CRR-E. Cortisol and glucose: reliable indicators of fish stress? 2009.Google Scholar
  10. Iwama GK, Vijayan MM, Forsyth RB, Ackerman PA. Heat shock proteins and physiological stress in fish. Am Zool. 1999;39(6):901–9.View ArticleGoogle Scholar
  11. Gracey AY. Interpreting physiological responses to environmental change through gene expression profiling. J Exp Biol. 2007;210(Pt 9):1584–92.View ArticlePubMedGoogle Scholar
  12. Teigen LE, Orczewska JI, McLaughlin J, O’Brien KM. Cold acclimation increases levels of some heat shock protein and sirtuin isoforms in threespine stickleback. Comp Biochem Physiol A Mol Integr Physiol. 2015;188:139–47.View ArticlePubMedGoogle Scholar
  13. Windisch HS, Frickenhaus S, John U, Knust R, Portner HO, Lucassen M. Stress response or beneficial temperature acclimation: transcriptomic signatures in Antarctic fish (Pachycara brachycephalum). Mol Ecol. 2014;23(14):3469–82.View ArticlePubMedGoogle Scholar
  14. Ju Z, Dunham RA, Liu Z. Differential gene expression in the brain of channel catfish (Ictalurus punctatus) in response to cold acclimation. Mol Genet Genomics. 2002;268(1):87–95.View ArticlePubMedGoogle Scholar
  15. Long Y, Song G, Yan J, He X, Li Q, Cui Z. Transcriptomic characterization of cold acclimation in larval zebrafish. BMC Genomics. 2013;14:612.View ArticlePubMedPubMed CentralGoogle Scholar
  16. Hori TS, Gamperl AK, Booman M, Nash GW, Rise ML. A moderate increase in ambient temperature modulates the Atlantic cod (Gadus morhua) spleen transcriptome response to intraperitoneal viral mimic injection. BMC Genomics. 2012;13:431.View ArticlePubMedPubMed CentralGoogle Scholar
  17. Hassinen M, Paajanen V, Haverinen J, Eronen H, Vornanen M. Cloning and expression of cardiac Kir2.1 and Kir2.2 channels in thermally acclimated rainbow trout. Am J Physiol Regul Integr Comp Physiol. 2007;292(6):R2328–39.View ArticlePubMedGoogle Scholar
  18. Deane EE, Woo NY. Cloning and characterization of the hsp70 multigene family from silver sea bream: modulated gene expression between warm and cold temperature acclimation. Biochem Biophys Res Commun. 2005;330(3):776–83.View ArticlePubMedGoogle Scholar
  19. Ewart KV, Richards RC, Driedzic WR. Cloning of glycerol-3-phosphate dehydrogenase cDNAs from two fish species and effect of temperature on enzyme expression in rainbow smelt (Osmerus mordax). Comp Biochem Physiol B Biochem Mol Biol. 2001;128(3):401–12.View ArticlePubMedGoogle Scholar
  20. Spence R, Gerlach G, Lawrence C, Smith C. The behaviour and ecology of the zebrafish, Danio rerio. Biol Rev Camb Philos Soc. 2008;83(1):13–34.View ArticlePubMedGoogle Scholar
  21. 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–73.View ArticlePubMedGoogle Scholar
  22. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136(2):215–33.View ArticlePubMedPubMed CentralGoogle Scholar
  23. Cross CE, Tolba MF, Rondelli CM, Xu M, Abdel-Rahman SZ. Oxidative stress alters miRNA and gene expression profiles in villous first trimester trophoblasts. Biomed Res Int. 2015;2015:257090.View ArticlePubMedPubMed CentralGoogle Scholar
  24. Ci D, Song Y, Tian M, Zhang D. Methylation of miRNA genes in the response to temperature stress in Populus simonii. Front Plant Sci. 2015;6:921.View ArticlePubMedPubMed CentralGoogle Scholar
  25. Chang TC, Wentzel EA, Kent OA, Ramachandran K, Mullendore M, Lee KH, Feldmann G, Yamakuchi M, Ferlito M, Lowenstein CJ, et al. Transactivation of miR-34a by p53 broadly influences gene expression and promotes apoptosis. Mol Cell. 2007;26(5):745–52.View ArticlePubMedPubMed CentralGoogle Scholar
  26. Leung AK, Sharp PA. microRNAs: a safeguard against turmoil? Cell. 2007;130(4):581–5.View ArticlePubMedGoogle Scholar
  27. Leung AK, Sharp PA. MicroRNA functions in stress responses. Mol Cell. 2010;40(2):205–15.View ArticlePubMedPubMed CentralGoogle Scholar
  28. 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–27.View ArticlePubMedPubMed CentralGoogle Scholar
  29. Yan B, Zhao LH, Guo JT, Zhao JL. miR-429 regulation of osmotic stress transcription factor 1 (OSTF1) in tilapia during osmotic stress. Biochem Biophys Res Commun. 2012;426(3):294–8.View ArticlePubMedGoogle Scholar
  30. Yang R, Dai Z, Chen S, Chen L. MicroRNA-mediated gene regulation plays a minor role in the transcriptomic plasticity of cold-acclimated zebrafish brain tissue. BMC Genomics. 2011;12:605.View ArticlePubMedPubMed CentralGoogle Scholar
  31. 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.View ArticlePubMedPubMed CentralGoogle Scholar
  32. Lukowski SW, Fish RJ, Martin-Levilain J, Gonelle-Gispert C, Buhler LH, Maechler P, Dermitzakis ET, Neerman-Arbez M. Integrated analysis of mRNA and miRNA expression in response to interleukin-6 in hepatocytes. Genomics. 2015;106(2):107–15.View ArticlePubMedGoogle Scholar
  33. Okamura H. Clock genes in cell clocks: roles, actions, and mysteries. J Biol Rhythms. 2004;19(5):388–99.View ArticlePubMedGoogle Scholar
  34. Yu W, Nomura M, Ikeda M. Interactivating feedback loops within the mammalian clock: BMAL1 is negatively autoregulated and upregulated by CRY1, CRY2, and PER2. Biochem Biophys Res Commun. 2002;290(3):933–41.View ArticlePubMedGoogle Scholar
  35. Rodriguez A, Griffiths-Jones S, Ashurst JL, Bradley A. Identification of mammalian microRNA host genes and transcription units. Genome Res. 2004;14(10A):1902–10.View ArticlePubMedPubMed CentralGoogle Scholar
  36. 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.View ArticlePubMedPubMed CentralGoogle Scholar
  37. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.View ArticleGoogle Scholar
  38. Long Y, Li L, Li Q, He X, Cui Z. Transcriptomic characterization of temperature stress responses in larval zebrafish. PLoS One. 2012;7(5):e37209.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Kanehisa M, Goto S, Furumichi M, Tanabe M, Hirakawa M. KEGG for representation and analysis of molecular networks involving diseases and drugs. Nucleic Acids Res. 2010;38(Database issue):D355–60.View ArticlePubMedGoogle Scholar
  40. Dunlap JC. Molecular bases for circadian clocks. Cell. 1999;96(2):271–90.View ArticlePubMedGoogle Scholar
  41. Honma S, Kawamoto T, Takagi Y, Fujimoto K, Sato F, Noshiro M, Kato Y, Honma K. Dec1 and Dec2 are regulators of the mammalian molecular clock. Nature. 2002;419(6909):841–4.View ArticlePubMedGoogle Scholar
  42. Laposky AD, Bass J, Kohsaka A, Turek FW. Sleep and circadian rhythms: key components in the regulation of energy metabolism. FEBS Lett. 2008;582(1):142–51.View ArticlePubMedGoogle Scholar
  43. Imaizumi T, Kay SA, Schroeder JI. Circadian rhythms. Daily watch on metabolism. Science. 2007;318(5857):1730–1.View ArticlePubMedPubMed CentralGoogle Scholar
  44. Wang Q, Tan X, Jiao S, You F, Zhang PJ. Analyzing cold tolerance mechanism in transgenic zebrafish (Danio rerio). PLoS One. 2014;9(7):e102492.View ArticlePubMedPubMed CentralGoogle Scholar
  45. Scott GR, Johnston IA. Temperature during embryonic development has persistent effects on thermal acclimation capacity in zebrafish. Proc Natl Acad Sci U S A. 2012;109(35):14247–52.View ArticlePubMedPubMed CentralGoogle Scholar
  46. Wang H, Zou Z, Wang S, Gong M. Global analysis of transcriptome responses and gene expression profiles to cold stress of Jatropha curcas L. PLoS One. 2013;8(12):e82817.View ArticlePubMedPubMed CentralGoogle Scholar
  47. Chou MY, Hsiao CD, Chen SC, Chen IW, Liu ST, Hwang PP. Effects of hypothermia on gene expression in zebrafish gills: upregulation in differentiation and function of ionocytes as compensatory responses. J Exp Biol. 2008;211(Pt 19):3077–84.View ArticlePubMedGoogle Scholar
  48. Simon-Oppermann C, Hammel HT, Simon E. Hypothalamic temperature and osmoregulation in the Pekin duck. Pflugers Arch. 1979;378(3):213–21.View ArticlePubMedGoogle Scholar
  49. Shin HS, Bargiello TA, Clark BT, Jackson FR, Young MW. An unusual coding sequence from a Drosophila clock gene is conserved in vertebrates. Nature. 1985;317(6036):445–8.View ArticlePubMedGoogle Scholar
  50. Ostberg O. Circadian rhythms of food intake and oral temperature in “morning” and “evening” groups of individuals. Ergonomics. 1973;16(2):203–9.View ArticlePubMedGoogle Scholar
  51. Lahiri K, Vallone D, Gondi SB, Santoriello C, Dickmeis T, Foulkes NS. Temperature regulates transcription in the zebrafish circadian clock. PLoS Biol. 2005;3(11):e351.View ArticlePubMedPubMed CentralGoogle Scholar
  52. Zani F, Breasson L, Becattini B, Vukolic A, Montani JP, Albrecht U, Provenzani A, Ripperger JA, Solinas G. PER2 promotes glucose storage to liver glycogen during feeding and acute fasting by inducing Gys2 PTG and G L expression. Mol Metab. 2013;2(3):292–305.View ArticlePubMedPubMed CentralGoogle Scholar
  53. Kennaway DJ, Varcoe TJ, Voultsios A, Boden MJ. Global loss of bmal1 expression alters adipose tissue hormones, gene expression and glucose metabolism. PLoS One. 2013;8(6):e65255.View ArticlePubMedPubMed CentralGoogle Scholar
  54. So AY, Bernal TU, Pillsbury ML, Yamamoto KR, Feldman BJ. Glucocorticoid regulation of the circadian clock modulates glucose homeostasis. Proc Natl Acad Sci U S A. 2009;106(41):17582–7.View ArticlePubMedPubMed CentralGoogle Scholar
  55. Wendelaar Bonga SE. The stress response in fish. Physiol Rev. 1997;77(3):591–625.PubMedGoogle Scholar
  56. Grimaldi B, Bellet MM, Katada S, Astarita G, Hirayama J, Amin RH, Granneman JG, Piomelli D, Leff T, Sassone-Corsi P. PER2 controls lipid metabolism by direct regulation of PPARgamma. Cell Metab. 2010;12(5):509–20.View ArticlePubMedPubMed CentralGoogle Scholar
  57. Kultz D. Molecular and evolutionary basis of the cellular stress response. Annu Rev Physiol. 2005;67:225–57.View ArticlePubMedGoogle Scholar

Copyright

© The Author(s). 2016

Advertisement