Global identification of the genetic networks and cis-regulatory elements of the cold response in zebrafish

Abstract

The transcriptional programs of ectothermic teleosts are directly influenced by water temperature. However, the cis- and trans-factors governing cold responses are not well characterized. We profiled transcriptional changes in eight zebrafish tissues exposed to mildly and severely cold temperatures using RNA-Seq. A total of 1943 differentially expressed genes (DEGs) were identified, from which 34 clusters representing distinct tissue and temperature response expression patterns were derived using the k-means fuzzy clustering algorithm. The promoter regions of the clustered DEGs that demonstrated strong co-regulation were analysed for enriched cis-regulatory elements with a motif discovery program, DREME. Seventeen motifs, ten known and seven novel, were identified, which covered 23% of the DEGs. Two motifs predicted to be the binding sites for the transcription factors Bcl6 and Jun, respectively, were chosen for experimental verification, and they demonstrated the expected cold-induced and cold-repressed patterns of gene regulation. Protein interaction modeling of the network components followed by experimental validation suggested that Jun physically interacts with Bcl6 and might be a hub factor that orchestrates the cold response in zebrafish. Thus, the methodology used and the regulatory networks uncovered in this study provide a foundation for exploring the mechanisms of cold adaptation in teleosts.

Figure 1.

Characterization of the cold-induced gene expression patterns in zebrafish tissues. (A) Schematic diagram showing the cooling time course and the sampling regimes. The fish were subjected to the following three temperatures: 28, 18 and 10°C. The red dots indicate the sampling time points. At each point, 20 fish were sampled, and RNA from eight tissues was isolated and pooled by tissue type for RNA-Seq. (B) The number of RNA-Seq reads obtained from the eight analysed tissues at the three temperatures. (C) The read distributions among the annotated zebrafish genomic features. (D) The tissue distributions of DEGs with a >2-fold RPKM difference for the three comparisons (28°C compared with 18, 18°C compared with 10 and 28°C compared with 10°C). The sizes of the circles are proportional to the numbers of genes that they contain. (E) MDS plot showing the diversity of the responses to the cold treatment among the zebrafish tissues. The entire set of expressed genes (22 457 genes) was used in analysis. The MDS algorithm placed each sample in a 2-dimensional space for visualization of the level of similarity of the tissue expression patterns. Each sample was assigned coordinates based on the expression levels of all genes in the sample, and these coordinates were used to construct a scatterplot. The three temperatures are indicated by different colours. Abbreviations: br for brain, gi for gill, he for heart, in for intestine, ki for kidney, li for liver, mu for muscle and sp for spleen.

Figure 2.

Comparison of the expression patterns of selected genes detected in RNA-Seq and RT-qPCR assays showing high correlation between the two methods. The log₂ ratios of the expression changes at the low temperatures relative to the normal temperature were calculated and plotted (the ratio was set to 0 for the normal condition). The red line depicts the RT-qPCR results. The expression levels of the selected genes were normalized against that of β-actin. The blue line shows the RNA-Seq results. The R² values (Pearson correlation coefficients) across the different temperatures are shown for each gene. The expression patterns of three genes—cirbp (A), atf4b2 (B) and id1 (C)—showed excellent agreement between the RNA-Seq and RT-qPCR assays. A boxplot of R² values for the 12 selected genes shows a high R² value (>0.7) across the tissues examined (D). Supplementary Figures S4 shows 9 additional genes that were also verified.

Figure 3.

Heatmap showing the enriched GO categories within the gene clusters. All clusters listed in Supplementary Table S3 were examined for overrepresented GO terms by comparisons with the entire expressed gene set. Enriched GO terms were identified for the 16 clusters shown at the top. The associated biological process GO terms identified in each cluster are plotted in the bottom heatmap panel. The genes annotated with each enriched GO term are shown to the right of the panel. For enriched GO terms associated with genes from two or more clusters, the cluster(s) to which the genes belong are denoted by Cluster IDs in brackets. The color scales depict the p values for the enrichment tests, and the grey cells indicate a P value of >0.05. The results of the GO enrichment test are shown in Supplementary Table S5.

Figure 4.

Identification of cold-responsive cis-regulatory motifs and distribution of the motif-containing genes among the enriched GO terms. (A) Strategies for identifying enriched motifs in the proximal promoter regions of the co-regulated genes. Genes were assigned to a cluster at two different membership cut-off values—0.12 and 0.4—based on the Pearson correlation of the gene's expression pattern with the centroid of the cluster. Conserved motifs that were within 1 kb upstream of their annotated 5′ UTRs were identified using the DREME algorithm. TOMTOM was used to annotate the detected motifs with JASPAR database. (B) Heatmap showing the percentages of genes in the gene set with certain motifs (shown at the top) associated with the enriched GO terms (shown at right) identified in Figure 3.

Figure 5.

Experimental verification of the regulatory roles of two motifs, AGMAACCA and SAGTCAA, in gene expression in response to cold stress. (A and B) The locations, sequences and orientations of the cis-motifs within the gene promoters cloned for transfection assays. The highly conserved motif sequences in each promoter are indicated by the gray boxes. (C) The scheme for constructing a plasmid containing the luciferase reporter gene driven by the native promoter (top) or the motif-replaced promoter (bottom). (D) Ratios of reporter mRNA induction or repression at 10°C relative to 28°C were analysed using the wild-type and mutant constructs by transfecting zebrafish ZF4 and grass carp CIK cell lines. Statistically significant differences in mRNA induction or repression between the wild-type and mutant promoters are indicated by asterisks (* indicates P < 0.05, and ** indicates P < 0.01, t-test). An empty vector, PGL4.10, was used as a control. WT: wild-type. (E) The experimental design for examining the cold-responsive functions of the selected motifs by transgenic study of zebrafish embryos. The structure of the promoter construct and the workflow are shown. (F) Green fluorescence images, showing EGFP expression under the control of the wild-type (top) and mutant (bottom) gstp1 promoters. (G) Relative expression levels of egfp in transgenic zebrafish embryos at cold (10°C) and normal (28°C) temperatures, as determined by RT-qPCR. The wild-type and mutated promoters of the same three genes (gstp1, ubc and apoba) as those used in the cell transfection experiments were examined. Fold changes were calculated by normalizing the mRNA levels at 10°C to those at 28°C. Student's t test was used to assess significant differences between the mutant/WT ratio and a value of 1. (H) Model of cis-regulation of cold-induced (red) and cold-repressed (blue) gene expression by Bcl6 and AP-1, respectively. Bcl6 positively regulates 39 genes from Cluster 0 via the AGMAACCA cis-motif, whereas AP-1 down-regulates the expression of 80 genes via the SAGTCAA motif. The known binding sites for Bcl6 and AP-1 are indicated at the top of each panel. The region that was predicted to contain the enriched motif by DREME is underlined.

Figure 6.

Regulatory networks, showing the protein and gene interactions among the motif-bearing genes and the binding TFs. (A) Circular nodes denote the classifications of the genes in the identified gene clusters, with red indicating the AGMAACCA motif in Cluster 0, blue depicting the SAGTCAA motif in Cluster 19, magenta denoting the TTCCAGBA motif in Cluster 1, green depicting the WCACCTGW motif in Cluster 7 and cyan indicating the CCAATCAG motif in Cluster 21. The quadrate nodes denote the binding TFs, with each colour indicating the TF–promoter relationship with the gene of the same colour. Protein–protein interaction information was retrieved from the literature using GeneMANIA. The types of protein–protein interactions are depicted by different types of lines, with the solid lines representing physical interactions, the parallel lines depicting co-localization, and the dots representing genetic interactions. According to the definitions in GeneMANIA, ‘physical interaction’ refers to an experimentally validated protein–protein interaction, and ‘co-localization’ refers to two proteins that co-localize in the cell. ‘Genetic interaction’ indicates that the effects of one gene are modified by another gene, as identified through mutational analysis. The number of connections associated with each TF is indicated in the quadrate nodes. Protein–protein interaction information was retrieved from the literature using GeneMANIA, and the network was constructed using Cytoscape software. (B) Validation of the interaction between Jun and Bcl6 by CoIP. ZF4 cells maintained at 28°C and 10°C were harvested for immunoprecipitation (IP) with a Jun or Bcl6 antibody. Input (top) and immunoprecipitates of Jun (middle panel) and Bcl6 (bottom panel) were analysed by immunoblotting (IB) with antibodies against Jun and Bcl6.