Cold temperature represses daily rhythms in the liver transcriptome of a stenothermal teleost under decreasing day length

In the strongly seasonal habitats in high latitudes, organisms anticipate seasonal changes in temperature by tracking daily and annual variation in photoperiod (Wood and Loudon, 2014). Photoperiod regulates the central circadian clock, an oscillating transcriptional feedback loop located in mammals in the suprachiasmatic nucleus of the hypothalamus, and in fish in the pineal gland, although in fish the peripheral clocks may be more independent (Idda et al., 2012; Lowrey and Takahashi, 2004). The central clock acts as a master regulator of the rhythms of peripheral tissues (Idda et al., 2012; Lowrey and Takahashi, 2004), but in addition to central clock entrainment, nutrient and redox status and other factors can also entrain the clock (Asher and Schibler, 2011; Yang et al., 2006), and in fish the peripheral rhythms can be directly light entrained (Whitmore et al., 2000). Across tissues, transcriptional regulation is a central feature of circadian rhythms (Akhtar et al., 2002; Doherty and Kay, 2010). However, the genes under circadian control are tissue specific (Panda et al., 2002).

Besides photoperiod, temperature can regulate circadian rhythms. In a temperate fish, the Atlantic salmon (Salmo salar), acclimation at a warm temperature can increase melatonin secretion during the dark phase (Porter et al., 2001), and in zebrafish (Danio rerio), temperature bouts and daily cycles can entrain the clock (Lahiri et al., 2005, 2014; López-Olmeda and Sánchez-Vazquez, 2009). However, temperature compensation within an ecologically relevant temperature range is important for an endogenous clock in order to maintain a ca. 24-h period length, particularly in ectotherms (Brown and Webb, 1948). However, the mechanism by which this is achieved and whether the clocks of cold-adapted species conform to this principle are not well known.

Deeper insight into the biological clocks of species at high latitudes is now highly warranted, as the temperature increase and extreme weather events associated with climate change can introduce a mismatch between an organism's circannual clock and the environmental conditions, which is likely to have negative impacts on wild populations (Jørgensen and Johnsen, 2014; Stevenson et al., 2015). Particularly in ectotherms, increasing temperature could also lead to changes in circadian rhythms within seasons. Adverse climatic effects are especially emphasized in the Arctic (Marshall et al., 2014), and projections of climate change over the next decades show that warming in the Arctic will likely be strongest during the autumn and winter (Walsh, 2008). As circadian variation in transcription can mask the responses to temperature for many genes (Podrabsky and Somero, 2004), accounting for daily variation in gene expression studies can also lead to major improvements in understanding the mechanisms affecting temperature tolerance. This has particular importance for research on stenothermal species that do not cope well at high temperatures and whose adaptive responses we need to understand better, such as the Arctic char (Anttila et al., 2015; Jeuthe et al., 2015; Penney et al., 2014).

In plants, cold temperatures have a direct input into the circadian machinery to repress processes that would otherwise be triggered by seasonal photoperiod changes, such as flowering (Mizuno et al., 2014). Arctic ectotherms would also likely experience a maladaptive response if only photoperiod was involved in signalling life-history events. Therefore, we tested whether cold temperature during a period of late summer photoperiod would result in a different pattern of transcript expression compared with warm temperature, particularly with respect to transcripts that have a daily rhythm. Consequently, our hypothesis was that temperature affects transcriptional rhythmicity in a stenothermal fish during late summer day lengths that coincide with the annual temperature maximum in the wild. To our knowledge, this is the first study evaluating temperature effects on overall transcriptional rhythms of a mostly Arctic species. Therefore, the results form a starting point for more detailed and ecologically relevant studies.

The population of Arctic char [Salvelinus alpinus (Linnaeus 1758)] used in the experiments originated from Lake Kuolimo, Finland (61°16′N, 27°32′E). The experiments were conducted at the Natural Resources Institute Finland (Luke) in Enonkoski, eastern Finland, in 2013. One-year-old Arctic char from a brood stock that has been maintained at the aquaculture facilities for three generations were reared under conditions of a natural light–dark rhythm and ambient temperature (approximately 8°C before the acclimations) using water fed from Lake Pahkajärvi (62°04′N, 28°33′E) prior to the experiment. In early July, warm acclimation was initiated by transferring fish to a tank that was heated to 15°C (mean 14.9±0.6°C) over 1 day, while another group of fish was maintained at a cold temperature (mean 7.7±0.2°C, from hereon 8°C). The test temperatures were roughly at the low and high end of the temperatures experienced and tolerated by the fish during summer; the cold water was the coldest available in the lake in the summer, and negative consequences for fish health from temperatures higher than 15°C had previously been found in the hatchery. Notably, the temperature range 8–15°C can occur simultaneously at different depths in lakes populated by Arctic char. There was no daily variation in water temperature in either group. The acclimations were conducted under natural light conditions, the length of light period being 19 h 20 min at the beginning and 17 h 15 min at the end of the acclimation period. No artificial lighting was provided for the fish; instead, all light entered through windows located on two sides of the room. The tanks were not exposed to direct sunlight. Fish were acclimated in identical, round, approximately 320 litre flow-through tanks. Oxygen levels were kept above 80% air saturation throughout the acclimations, and temperatures of the tanks were monitored and adjusted daily. Mean O₂ levels during acclimation were 10.2±2.3 mg l⁻¹ in the warm treatment and 12.5±2.3 mg l⁻¹ in the cold treatment. Thus, we minimized random variation between the two tanks using controlled conditions (feeding, oxygen level, density and size of fish, uniform photoperiod and light intensity, removal of uneaten food and faeces). Fish were fed with commercial fish pellets (Raisio Group; www.raisiogroup.com) ad libitum using automatic feeders and fasted for 24 h before sampling.

After the 4-week acclimation, fish from each temperature were collected at three time points during the day – at 12:30 h (middle of the light period), 20:30 h (approximately 1 h before sunset) and 01:30 h (middle of the dark period), hereafter referred to as day, evening and night, respectively – euthanized with 200 ppm sodium bicarbonate-buffered tricaine methanesulfonate (MS-222), then measured and weighed. Blood was collected in heparinized capillary tubes for haematocrit measurements. Liver tissue was removed and flash-frozen in liquid nitrogen, after which the samples were stored at −80°C until analysis. The fish were sexed during the dissection by inspecting the gonads. The handling time for each fish from euthanasia to removal of tissue was <3 min.

All procedures were approved by the Finnish Animal Experiment Board (ESAVI/4068/04.10.07/2013).

Three males for each time point–temperature combination with roughly uniform size [body mass, 15°C=22.2±8.4 g (mean±s.d.), 8°C=24.2±10.4 g; fork length at 15 and 8°C=13.7±1.8 cm] were selected for RNA extraction and sequencing (total N=18). These fish also had uniform haematocrit values (47±4% and 46±5% at 8 and 15°C, respectively). RNA was extracted from approximately 10 mg liver tissue at the University of Turku. Extractions were performed using Tri Reagent (Molecular Research Center, Cincinnati, OH, USA) following the manufacturer's protocol. In short, tissue was homogenized using TissueLyser (Qiagen, Hilden, Germany) for 90 s at 30 Hz. Eluted RNA was treated with DNase I and purified with a re-extraction using Tri Reagent. RNA was quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and its integrity was confirmed using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) (average RNA integrity number 9.95). Strand-specific cDNA library preparation and sequencing were conducted at the Beijing Genomics Institute (BGI Hong Kong) using the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA) and an Illumina HiSeq2000 instrument. The reads were trimmed from adapter sequences at the sequencing facility. We checked the sequencing data for quality with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and trimmed reads using Trimmomatic (Bolger et al., 2014) (details in the Appendix).

The whole genome sequence of Arctic char is not yet available; thus, a de novo transcriptome assembly was generated from quality-trimmed paired-end (PE) reads using Trinity v.2.0.6 (Grabherr et al., 2011) (details in the Appendix). The trimmed PE reads from all samples were aligned against the assembly using Bowtie 2 (Langmead and Salzberg, 2012), ignoring discordant and mixed alignments and allowing a maximum of 40 multiple alignments per read. Read counts per cluster were obtained from generated alignment files using Corset (Davidson and Oshlack, 2014). Clusters that were expressed with >1 counts per million (cpm) reads in at least three samples across all 18 samples were included in the analysis, and are hereafter referred to as genes, although some of the assembled sequences may not correspond to actual mRNAs, but instead to long non-coding or pre-mRNAs. Reads for gene isoforms and paralogues of high similarity may be counted under the same clusters, but this issue remains challenging to circumvent in a de novo assembly analysis. The assembly was filtered to include only the transcripts of expressed genes to generate a final assembly available for download.

Differential expression (DE) analysis was conducted using the limma package v.3.28.21 (Ritchie et al., 2015; Smyth, 2005) in RStudio v. 0.99.903 (R v. 3.3.1, R Foundation for Statistical Computing, Vienna, Austria). The data were normalized for library sizes using the edgeR package TMM (trimmed mean of M-values) normalization (Robinson and Oshlack, 2010). The voom method was applied to obtain precision weights for the mean–variance relationships of expression levels (Law et al., 2014) that were included in the limma empirical Bayes analysis pipeline (Ritchie et al., 2015; Smyth, 2005). We built four contrast matrices for DE analysis: (1) the average temperature effect across time points, (2) pairwise differences between time points within temperatures, (3) the interaction effect of time point and temperature and (4) differences between temperatures within time points. The average effects of time points were not included in the analysis owing to the strong temperature dependence of the effect. The P-values were adjusted for false discovery rate (FDR) using the Benjamini–Hochberg method, and the adjusted P-value 0.05 was used as the significance threshold, calculated using the limma decideTest method ‘global’ across contrasts in matrices 2 and 4 to compare the strengths of different contrasts.

To evaluate the proportion of total variation between samples explained by time or temperature, the voom-transformed and normalized read counts were further used in partial redundancy analysis (RDA) using the package vegan v. 2.4-3 (https://CRAN.R-project.org/package=vegan). Variation between samples was explained by: (1) time, temperature and their interaction, (2) time after conditioning for the effect of temperature and (3) temperature after conditioning the effect of time. Clustering of samples was visualized using scaling based on genes.

The normalized log₂ cpm data generated by voom for selected target genes were visualized using the packages plyr, ggplot2 and gplots in R (Warnes et al., 2016; Wickham, 2009, 2011). R codes used in the analysis and visualization of read counts are provided in Script 1.

Open reading frame (ORF) peptide sequences were obtained for transcripts in the final assembly using TransDecoder software (https://github.com/TransDecoder/TransDecoder/wiki). The predicted ORFs were annotated with three databases using the Basic Local Alignment Search Tool for proteins (BLASTp v.2.2.31). All of the predicted ORF sequences were annotated with predicted zebrafish (downloaded 24 October 2015 from Ensembl) and salmon (NCBI Salmo salar Annotation Release 100) proteins using a reciprocal best hits approach and an e-value cut-off of 1×10⁻⁵. Additionally, the ORFs were annotated with the NCBI non-redundant protein database (NCBI nr, downloaded 25 November 2015) with an e-value cut-off of 1×10⁻⁵ and when the query sequence matched the target sequence at >50% protein length. The annotations were prioritized with the order zebrafish>salmon>NCBI nr. When available, gene descriptions were retrieved from Ensembl using biomaRt in R. Gene symbols were retrieved for zebrafish Ensembl IDs, salmon Refseq IDs and NCBI gene names using the Biological DataBase Network (https://biodbnet-abcc.ncifcrf.gov). Annotations were retrieved for 9491, 4037 and 4117 genes with zebrafish peptides, Atlantic salmon predicted peptides and the NCBI nr-database, respectively. In total, 20,394 out of 44,784 ORFs in the assembly (45.5%) were annotated with 18,013 unique protein IDs.

Because only a limited proportion of the annotations retrieved for the transcripts could be transformed into human IDs for further analyses, we also annotated all ORF peptide sequences with human peptide sequences using BLASTp v.2.4.0 with an e-value threshold of 10⁻⁵. After identifying human orthologues, we supplemented the annotation with the previously obtained gene symbols for genes that were missing an annotation. In total, 18,232 genes were annotated using this approach with 9577 unique gene symbols.

In order to investigate patterns of co-expression among transcripts, we also analyzed the normalized, voom-transformed data using weighted correlation network analysis (also known as weighted gene co-expression network analysis; WGCNA) (Langfelder and Horvath, 2008). We first removed transcripts with low variation across treatments (cpm variance <1) and ran the function blockwiseModules to identify potentially co-regulated genes. We created a signed network using a soft thresholding power of 20, a minimum module size of 30 transcripts and a merge cut height of 0.25. As it is likely that modules of genes are regulated by common transcription factors, we examined the annotated members of the 10 modules whose eigengenes were significantly correlated with temperature treatment or time using the upstream regulator analysis in Ingenuity Pathway Analysis (IPA) software (Qiagen). This analysis uses a Fisher's exact test to determine whether the input list has significantly more members that are known targets of transcription factors than a random list of genes.

We obtained on average 40.1±6.5 million (mean±s.d.) raw PE reads per sample (minimum 27.3 million reads, maximum 60.6 million reads) for the DE analysis. A total of 99% of reads passed quality trimming and were used in the assembly and read alignment. The assembly contained approximately 210,000 transcripts prior to and 88,487 transcripts after filtering based on expression values. The transcripts originated from approximately 42,000 genes. The alignment rate was on average 23.0±6.0% for unique and 60.3± 1.2% for multiple-aligned reads. Inspection of sam-files generated by Bowtie 2 revealed that multiple-aligned reads were generally aligning to different isoforms of the same Trinity genes, thus having a small effect on the gene-level expression analysis. As determining isoform-specific expression accurately was not possible because of the short read length and presence of multiple paralogous genes (as a result of a salmonid-specific whole genome duplication event; e.g. Macqueen and Johnston, 2014), we focused on the gene-level differences in the DE analysis. However, many conserved gene paralogues of salmonids were found as potentially separate transcripts or genes in the assembly (Table S1). BUSCO analysis using conserved, near-universal single-copy orthologues (Simão et al., 2015) indicated that the assembly had 61% completeness, 4.7% fragmentation and 33% missing orthologues of these genes. The transrate (Smith-Unna et al., 2016) score of the assembly was increased by expression-level filtering from 0.18 to 0.24.

Clustering of samples using partial RDA showed the samples from different temperatures could be clearly separated by their pairwise distances (Fig. 1A), and that much larger distances separated samples from different time points at 15°C than at 8°C after conditioning for the effect of temperature (Fig. 1A). Likewise, the samples were clearly separated between temperatures after conditioning for the effect of time, with day and evening samples clustering together within time points (Fig. 1C). In the DE analysis, we found significant pairwise differences between time points in 2080 genes at 15°C (Table S2), representing ca. 5% of the whole transcriptome, whereas 41 genes showed temporal variation at 8°C (Table S3), only five of which were annotated. Nine genes showed temporal variation at both temperatures (two annotated: filamin A interacting protein 1-like and zinc finger protein 2-like). At 15°C, the differences in gene expression between time points were largest between the day and evening, and the direction of change from day to evening or night was, for the most part, an increase in expression (Table 1; Table S2, Fig. S1). In addition, we found 2925 genes showing a significant time–temperature interaction effect (Table S4).

We detected differential expression in two orthologues of circadian locomotor output cycles kaput 1 (clock1). A comparison of the cDNA sequences of clock genes across vertebrates revealed one of the clock transcripts to be more similar to clock1a (cluster 92957.0), while the other was more similar to clock1b (cluster 78073.3); it is referred to as clock1b from here on. Clock1a (Fig. S2A), as well the orthologues of aryl hydrocarbon receptor nuclear translocator-like protein 2 (arntl2, also known as bmal/cycle2) and cryptochrome 1aa (cry1aa) were more highly expressed at 8°C than at 15°C. Clock1b was expressed less, and arntl1 expressed more, during the night and evening than during the day at 15°C (Fig. S2B,C). We additionally compared the differentially expressed genes with those from the super pathway ‘BMAL1-CLOCK, NPAS2 activates circadian gene expression’ in Path Cards, a pathway unification database (http://pathcards.genecards.org; Belinky et al., 2015), which lists 78 genes that includes the core clock components and their interaction partners. In total, we found that 36 members of this pathway were significantly differentially expressed between time points, temperatures or an interaction between the two (Table 2).

We contrasted the average effect of temperature across all time points in the DE analysis, which resulted in 4099 genes with higher expression in the warm than in the cold treatment, and 4633 genes with higher expression in the cold than in the warm treatment with FDR<0.05 (in total, 20% of the transcriptome; Fig. S3, Table S5). However, when comparing temperatures in each time point independently, the number of differentially expressed genes varied greatly (Fig. 2). Only 176 genes that were differentially expressed between temperatures were shared across all time points.

WGCNA clusters genes with similar patterns of expression across samples to create modules of genes that are likely co-expressed. Because this method uses hierarchical clustering of expression values to group genes into modules, the connectivity of the genes in the modules could reflect the response to time, temperature or their interaction. After the modules were created, the correlation of the module eigengenes with time and temperature was calculated to examine the strength of the correlation of the module with a given trait. There were two modules that were significantly correlated with time and eight modules that were significantly correlated with temperature (Table 3), although because the clustering is independent of the traits, modules could also reflect the interaction between traits. Details of the genes belonging to the modules are provided in Table S6.

When the expression of the genes within the modules is visualized using heat maps, the yellow module shows a very clear interaction between time and temperature, with the genes from the cold-acclimated individuals being uniformly expressed across all time points (with one outlier) and the warm-acclimated individuals showing a clear time response (Fig. 3). The green, turquoise and blue modules also show this interaction to a large extent, but there is slightly more variation among individuals (Figs S4–S6). The brown module was mainly clustering genes that are responding to the temperature difference but had no daily rhythm (Fig. S7). Significantly expressed genes from the DE analysis also overlapped with many genes in the modules correlated with time or temperature, whereas only seven significant DE genes belonged to modules not significantly correlated with those traits (Table 4).

In the turquoise module, which was significantly correlated with temperature but contained many genes that were significantly affected by time and the time–temperature interaction in the DE analysis, HNF4α was the top-scoring potential upstream regulator (Table S7). The expression of hnf4α was affected by the time–temperature interaction (Fig. S2D). There were also 19 of the 78 genes from the ‘BMAL1-CLOCK, NPAS2 activates circadian gene expression’ super pathway that were identified as potential upstream regulators of the genes in this module, 12 of which were also significant in the DE analysis results (Table 2). In comparison, only five genes from this super pathway were found among the upstream regulators of the yellow module, which was primarily correlated with time, highlighting the temperature-dependent effects of circadian regulators.

The main finding in our study was the profound thermal effect on daily rhythms in transcription, which was observed during a season of decreasing day length. In invertebrates, very low temperatures can stop the circadian clock from cycling (Brown and Webb, 1948; Hastings and Sweeney, 1957), and in ectothermic vertebrates, low temperature can dampen circadian rhythms in melatonin secretion or clock gene expression (Rensing and Ruoff, 2002; Vallone et al., 2007), but this study is the first to report a genome-wide decrease in transcriptional rhythms that is due to the combination of decreasing day length and cold temperature in fish. However, in Arctic char, melatonin rhythm persists throughout the year, suggesting that the pineal gland functions also at cold temperatures (Strand et al., 2008). A potentially similar phenomenon, where transcriptional rhythms of circadian clock genes can be paused by cold treatment, is known in plants (Martino-Catt and Ort, 1992). Therefore, we now have an indication that this response, where a decrease in temperature reduces or stops the cycling of specific transcripts, may be common for many different taxa. Having said the above, we note that a shift in the phase of the daily rhythm instead of decrease in amplitude may explain the results for at least some of the genes. However, this alternative is, in our opinion, less likely in the context of previous studies in other organisms (see above), and we would expect to see more genes showing significant fold-changes between time points at 8°C if the phase shift explained the majority of the findings. Our results also highlight that daily variation should be taken into consideration in future studies that apply genomics techniques to uncover adaptive responses, supporting previous observations (Tauber and Kyriacou, 2008).

Large daily variation in transcription and other processes can be mediated by endocrine signalling, as outlined in Cowan et al. (2017) for fishes. For instance, melatonin may affect the temperature-dependent transcriptional rhythms that we observed, for it plays a role in circadian regulation and warm temperature can increase night time melatonin secretion, as demonstrated in Atlantic salmon and white sucker (Catostomus commersoni) (Porter et al., 2001; Zachmann et al., 1992). Further, activity of the melatonin synthesis-limiting enzyme of teleosts, arylalkylamine N-acetyltransferase 2 (AANAT2), is correlated with the preferred temperature in a range of fishes, suggesting that there is an adaptive relationship between temperature and melatonin synthesis (Falcón et al., 2009). Other metabolic and endocrine responses may also play a role and can have joint effects on the temperature-dependent rhythms we observed. However, it is also possible that the effects on liver transcriptional rhythms are independent of the pineal gland, because melatonin secretion in Arctic char seems to occur also in cold temperatures (Strand et al., 2008).

In ectotherms, core clock gene expression shows changing amplitudes and peaks, but a conserved period length of 24 h, across temperatures (Kidd et al., 2015; Lahiri et al., 2005; Rensing and Ruoff, 2002), or a universal increase, accompanied by dampening of rhythms (Vallone et al., 2007). Likewise, cold temperature can dampen transcriptional rhythms in clock genes in Arabidopsis thaliana, a plant model species (Bieniawska et al., 2008). Our results showed that transcription of none of the core clock genes showed daily variation at 8°C, whereas two of them were rhythmic at 15°C (clock1b, arntl1) and five of them were more highly expressed at 8°C than at 15°C (clock1a, arntl1, arntl2, cry1aa, per1). Thus, the different clock genes may have divergent roles in integrating the thermal responses with circadian rhythms. Clock, arntl and cry are conserved members of the vertebrate circadian clock feedback loop, which oscillates in peripheral tissues and in the central clock with a 24-h rhythm. Cry transcription, along with the clock gene period and hundreds of other genes, is regulated by CLOCK and ARNTL binding on E-box elements, and the accumulation of PERIOD and CRY proteins gradually inhibits the DNA binding of CLOCK and ARNTL (Lowrey and Takahashi, 2004). In teleosts, different paralogues of circadian clock genes can have diverged roles in the feedback loop, as shown in zebrafish (Vatine et al., 2011). To our knowledge, there is no previous research on clock gene expression in Arctic char. Although it is unclear how temperature acts to repress the core clock components and how their protein levels are affected, it is likely that the lack or decrease of daily cycling of many downstream genes is the result of repressed rhythmicity of these core regulators, as studies in mice have identified thousands of clock-controlled genes in the liver tissue (e.g. Yoshitane et al., 2014; Rey et al., 2011). In future studies, comparing CLOCK protein accumulation at different acclimation temperatures and across populations with different thermal histories would be useful for understanding how it can trigger pathways related to shifts in life-history traits. CLOCK may be associated with seasonal regulation; the gene was found in a quantitative trait locus region associated with reproduction and growth in salmonids (Leder et al., 2006; Paibomesai et al., 2010) and clock alleles have a latitudinal gradient in the Chinook salmon (Oncorhynchus tshawytscha) and the blue tit (Cyanistes caeruleus) (Johnsen et al., 2007; O'Malley and Banks, 2008).

Owing to the small sample size for each time point in this study, the statistical analysis likely somewhat suffered from low power to detect differences between time points using the pairwise contrasts. However, the partial RDA and WGCNA results support the conclusions of the DE analysis. All three analyses indicated clear interaction between temperature and daily rhythms in regulating gene expression, with more pronounced changes at 15°C compared with 8°C. These changes could be mediated by transcriptional upstream regulators. Upstream analysis of the WGCNA module genes notably showed that hepatocyte nuclear factor 4α (HNF4α) was among the most significant upstream regulators of most of the modules tested in IPA. HNF4α binds to the promoters of almost half of the genes occupied by RNA polymerase II in human hepatocytes (Odom et al., 2004). It is a likely candidate for mediating the signalling between thermal (metabolic) regulation and the circadian clock, as disrupted light–dark cycles alter the transcription of HNF4α, as well as other hepatic gluconeogenic regulatory genes in mice (Oishi and Itoh, 2013). In our data, hnf4α expression was time dependent in contrasting directions at different temperatures (Fig. S2D).

Our gene expression data can provide new insights into the molecular mechanisms linking circadian regulation of transcription to organismal responses such as hypoxia and temperature tolerance in fish. Notably, we have previously shown that Arctic char acclimated to temperatures similar to those in the present study have significant differences in thermal tolerance, indicating physiological acclimation to the temperature treatments (Anttila et al., 2015). Further, understanding how temperature affects the patterns of daily variation in transcriptome expression can advance the development of chronotherapeutic approaches to treating fish diseases, as the toxicity of therapeutic agents, such as H₂O₂, can be rhythmic (Vera and Migaud, 2016). For instance, based on transcriptional patterns, daily variation in drug hepatotoxicity is expected to be low at relatively cold temperatures, but could vary throughout the day at warmer temperatures.

We found that daily rhythms in gene expression were dampened at 8°C but strong at 15°C in Arctic char, which suggests that daily transcriptional rhythms are strongly affected by temperature in this species, and leads to new hypotheses on the role of temperature in the regulation of biological clocks in fish. Firstly, the temperature compensation feature of circadian rhythms may have been under stronger selection in eurythermal than in stenothermal species, which generally experience small variation in temperature. Experiments on both steno- and eurythermal species under constant conditions at different temperatures could be used to address this. Secondly, temperature-dependent changes in the clock machinery may be induced by or enhance the negative effects of cold temperature on metabolic rate, because the circadian system interacts directly with metabolic regulators. Thirdly, circadian clock pathways could partially mediate the adjustment of life-history transitions by temperature, as they have previously been linked to the regulation of, for instance, smoltification and spawning in salmonids. The potentially diverged functional roles of circadian clock genes are of high interest in this respect.

Arctic char is known for its high phenotypic plasticity (Saether et al., 2015), and our study is focused on one sub-Arctic land-locked population, which may show divergent responses from Arctic populations. However, based on similar responses in expression observed in the clock genes in other fishes as well as in a phylogenetically very distant group, plants, it is likely that this response is conserved. In general, high phenotypic plasticity can contribute significantly to local adaptation, thereby promoting species' persistence in heterogeneous environments (Aubin-Horth and Renn, 2009; Dayan et al., 2015; Morris et al., 2014). In future studies, its role in the biological clocks of fish in relation to thermal adaptation should also be addressed.

The TruSeq RNA Sample Prep Kit v2 (Illumina) was used to construct libraries from 200 ng total RNA samples enriched for mRNA using oligo-dT beads. After RNA fragmentation, first-strand cDNA was generated by First Strand Master Mix and Super Script II (Invitrogen) reverse transcription (reaction conditions: 25°C for 10 min; 42°C for 50 min; 70°C for 15 min). The products were purified with Agencourt RNAClean XP Beads (Agencourt, Beckman Coulter, Brea, CA, USA). Second-strand cDNA synthesis was performed with Second Strand Master Mix and dNTP mix at 16°C for 1 h. The cDNA was combined with End Repair Mix, incubated at 30°C for 30 min and purified with Ampure XP Beads (Agencourt). PolyA-tails were added (A-Tailing Mix) and samples were incubated at 37°C for 30 min. RNA Index Adapter and Ligation Mix were added to the adenylated 3′ ends DNA and the samples were incubated at 30°C for 10 min. The end-repaired DNA was purified with Ampure XP Beads. To perform strand-specific sequencing, the second strand of cDNA was digested using the uracil-N-glycosylase enzyme and incubation at 37°C for 10 min, after which the samples were purified with Ampure XP Beads. Several rounds of PCR amplification with PCR Primer Cocktail and PCR Master Mix were performed to enrich the cDNA fragments. The PCR products were purified with Ampure XP Beads. The libraries were assessed for quality and quantity using two methods: checking the distribution of the fragments size using an Agilent 2100 bioanalyzer (Agilent DNA 1000 Reagents) and quantifying the library using real-time quantitative PCR (qPCR) with a TaqMan Probe. The quality-checked libraries were amplified on cBot to generate the cluster on the flow cell (TruSeq PE Cluster Kit V3-cBot-HS, Illumina). The amplified flow cell was sequenced using paired-end mode on the HiSeq 2000 System (TruSeq SBS KIT-HS V3, Illumina) using a 101 bp read length. Prior to sequencing, all samples were pooled and distributed across four lanes. Additional sequencing for the assembly construction was performed from the same samples pooled across two lanes.

The assembly and downstream analyses were conducted using computing resources from the IT Centre for Science, Espoo, Finland (www.csc.fi). Owing to challenges posed by accurately determining transcripts in a pseudotetraploid species, the transcriptome assembly was built in three steps. At first, the data from all individuals were combined and assembled using Trinity software v. 2.0.6 (Grabherr et al., 2011). The samples were aligned back to the transcriptome using Bowtie 2 v.2.2.5 (Langmead and Salzberg, 2012), allowing for a maximum of 40 multiple alignments when no unique alignment could be determined. Read counts per gene were analyzed with Corset software v.1.03, which clusters transcripts based on sequence similarity and expression levels (Davidson and Oshlack, 2014). Read count data were used to determine which samples covered the most variation in expression by determining the distances between samples based on all transcripts using the DeSeq2 package (Love et al., 2014) in R. The greatest distance was between several samples from 8°C and those from the evening time point at 15°C. Thus, the sample with the highest number of reads from both of these groups was chosen for the generation of the next assembly in order to reduce redundancy caused by using a large number of individuals. The second assembly was constructed using Trinity on trimmed reads obtained from two samples (in total, approximately 170 M read pairs) using a minimum kmer coverage of five reads and a minimum kmer length of 30 bp.

To improve the assembly, the samples used to generate the second assembly were aligned to the generated transcripts generated using Bowtie 2, and reads that were not concordantly aligned were used to generate a third assembly using Trinity. Finally, a BLAST search between the second and third assemblies was used to remove the redundant transcripts, and the transcripts from third assembly that had less than 98% similarity to those present in the second assembly and were longer than 100 bp were appended to the second assembly. Assembly quality was determined by the number of transcripts identified, rate of unique/multiple mapping reads, number of hits to the Uniprot Swissprot database using BLASTx, and the number of full-length hits from predicted open reading frames (ORFs) to the Uniprot Swissprot database using BLASTp and BUSCO, which estimates the completeness using universal single copy orthologues with the dataset for vertebrates (Simão et al., 2015). Additionally, the cDNA sequences from a set of 20 conserved paralogous genes in salmonids were compared with the final assembly using BLASTn.

leading:5 trailing:5 slidingwindow:4:15 minlen:36

–seqType fq –max_memory 96G –left all_25_pr.fq,all_56_pr.fq –right all_25_pf.fq,all_56_pf.fq –SS_lib_type RF –normalize_reads –min_kmer_cov 5 –kmer_size 30

–seqType fq –max_memory 96G –left 25_unmapped.1.fq,56_unmapped.1.fq\ –right 25_unmapped.2.fq,56_unmapped.2.fq –SS_lib_type RF –min_kmer_cov 2 –kmer_size 30

–phred64 –end-to-end -k 40 -D 20 –fr –no-mixed –no-discordant –score-min L,-0.1,-0.1 –no-unal

We appreciate the computing resources provided by the CSC-IT Center for Science, Espoo, Finland. We are grateful to Pasi Arkko for help in establishing the acclimations, to K. Field, B. J. G. Sutherland and F.-O. Gagnon-Hébert for discussions on the data analysis, and to D. Macqueen for providing the sequences of salmonid paralogues. We thank the anonymous reviewers for their constructive comments.