MicroRNAs regulate gene plasticity during cold shock in zebrafish larvae

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: Clock genes; Cold stress; Deep sequencing; Embryonic development; Zebrafish; miRNA.

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

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

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 (R² = 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

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

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

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

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)