Transcriptome Analysis of Hypothalamic Gene Expression during Daily Torpor in Djungarian Hamsters (Phodopus sungorus)

Animals living at high or temperate latitudes are challenged by extensive changes in environmental conditions over seasons. Djungarian hamsters (Phodopus sungorus) are able to cope with extremely cold ambient temperatures and food scarcity in winter by expressing spontaneous daily torpor. Daily torpor is a circadian controlled voluntary reduction of metabolism that can reduce energy expenditure by up to 65% when used frequently. In the past decades it has become more and more apparent, that the hypothalamus is likely to play a key role in regulating induction and maintenance of daily torpor, but the molecular signals, which lead to the initiation of daily torpor, are still unknown. Here we present the first transcriptomic study of hypothalamic gene expression patterns in Djungarian hamsters during torpor entrance. Based on Illumina sequencing we were able to identify a total number of 284 differentially expressed genes, whereby 181 genes were up- and 103 genes down regulated during torpor entrance. The 20 most up regulated group contained eight genes coding for structure proteins, including five collagen genes, dnha2 and myo15a, as well as the procoagulation factor vwf. In a proximate approach we investigated these genes by quantitative real-time PCR (qPCR) analysis over the circadian cycle in torpid and normothermic animals at times of torpor entrance, mid torpor, arousal and post-torpor. These qPCR data confirmed up regulation of dnah2, myo15a, and vwf during torpor entrance, but a decreased mRNA level for all other investigated time points. This suggests that gene expression of structure genes as well as the procoagulation factor are specifically initiated during the early state of torpor and provides evidence for protective molecular adaptions in the hypothalamus of Djungarian hamsters including changes in structure, transport of biomolecules and coagulation.


INTRODUCTION
Metabolic depression (torpor) is a commonly used strategy of mammals to survive winter. A reduction in energy expenditure as well as energy requirements is necessary to compensate for harsh environmental conditions during winter season when ambient temperature (T a ) drops and food availability is reduced (Jastroch et al., 2016).
The Djungarian hamster (also known as Siberian hamster, P. sungorus) has evolved a number of physiological and morphological adaptations (e.g., voluntary reduction of body weight, molt to dense white winter coat, gonadal regression, torpor) to seasonally reduce energy requirements (Figala et al., 1973;Scherbarth and Steinlechner, 2010). In Djungarian hamsters, winter adaptations are driven by photoperiod and can easily be induced by changes of the artificial light-dark cycle at moderate T a in the laboratory (Steinlechner and Heldmaier, 1982;Vitale et al., 1985). The most effective adaptive trait is the expression of daily torpor that spontaneously occurs after 10-12 weeks in short days once all other physiological adaptations are completed and the corresponding hormonal systems are in winter state (reduced levels of prolactin, testosterone and leptin) (Cubuk et al., 2016). Daily torpor is initiated by an active depression of metabolic rate (25% below the level of resting metabolic rate), accompanied by reduced heart rate and ventilation as well as a decrease in body temperature (T b ) to > 15 • C and reduced physical locomotor activity (Heldmaier and Ruf, 1992;Heldmaier et al., 2004). Torpor bouts are usually timed into the light phase of the light-dark cycle and have been shown to be under circadian control. The average duration of a torpor episode is 6 h and it is terminated by a spontaneous arousal prior to the hamsters' naturally active phase . The incidence of daily torpor is highly variable between individuals as well as in the same animal (1-7 torpor bouts per week) and can save up to 65% of total energy requirements, when torpor is used on a daily basis (Heldmaier, 1981;Kirsch et al., 1991;Ruf et al., 1991).
Spontaneous daily torpor is dependent on signaling of various hormonal systems changing with seasons and morphology, nutritional state as well as circadian timing mechanisms, hence, the hypothalamus is the brain area most likely involved in its regulation. Manipulations of prolactin levels lead to reduced torpor incidence and when testosterone, leptin or T3 are supplemented peripherally torpor is almost completely blocked Ruby et al., 1993;Freeman et al., 2004;Bank et al., 2015). It has already been shown, that lesion of various hypothalamic nuclei (suprachiasmatic nuclei, arcuate nucleus, paraventricular nucleus) alters torpor behavior. Moreover, the pharmacological activation of NPY mechanisms in arcuate nucleus induces torpor like hypothermia and hypothalamic application of T3 is able to block the expression of torpor (Ruby and Zucker, 1992;Ruby, 1995;Paul et al., 2005;Pelz and Dark, 2007;Dark and Pelz, 2008;Pelz et al., 2008). However, although torpor physiology has been extensively studied in this species, the regulatory systems in the brain ultimately initiating entrance into torpor on some days but not on others are entirely unknown.
Here we carried out a next generation sequencing (NGS) study to impartially screen for potential candidate genes playing a role in molecular hypothalamic torpor induction mechanisms. NGS allows the investigation of all transcripts of a genome by counting the number of mRNA sequencing reads of a specific tissue. To date, only few transcriptomic studies are available investigating gene expression patterns in the 13-lined ground squirrel (Ictidomys tridecemlieatus) during hibernation in various tissues, like cerebral cortex, hypothalamus, heart, skeletal muscle, brown adipose tissue and white adipose tissue (Hampton et al., 2011Schwartz et al., 2013Schwartz et al., , 2015Grabek et al., 2015). Except for one study investigating brown adipose tissue during entrance into torpor (Grabek et al., 2015), these studies were focused on time points before animals enter hibernation season, while being in deep hibernation or during the interbout arousals. The Djungarian hamster is an excellent animal model to investigate gene expression patterns during torpor entrance, because torpor is precisely timed into the circadian cycle and allows precise sampling with timed controls that are winter adapted but do not show torpor on that particular day. Moreover, substantial knowledge exists about hypothalamic mechanisms regulating seasonal adaptations in body weight and reproduction in this species (Ebling and Barrett, 2008).
Here we present a summary of differentially expressed genes during torpor entrance in the hypothalamus of P. sungorus. Moreover, we provide information about circadian regulation of mRNA expression patterns for selected candidate genes by relative gene expression analysis in torpid and normothermic hamsters.

Animals
All experiments and procedures were conducted in accordance with the German Animal Welfare Law and approved by the local animal welfare authorities (No. 114_14, Hamburg, Germany). All animals originated from our own breeding colony at the Institute of Zoology University of Hamburg. Djungarian hamsters (P. sungorus) were bred and raised under artificial long photoperiod (16:8-h light:dark cycle) at 21 • C ± 1 • C T a . The animals were individually housed in plastic cages (Macrolon Type III). Before and during the experiments, hamsters were fed a hamster breeding diet (Altromin 7014, Germany) ad libitum and had free access to drinking water. For the experiments, 3-4 months old Djungarian hamsters were transferred to short day conditions (8:16-h light:dark cycle) at constant T a of 18 • C ± 1 • C. After 12 weeks in short days they were implanted i.p. with DSItransmitters (Model TA-F10, St. Paul, MN, USA) under 1.5-2% isoflurane anesthesia and carprofen (5 mg/kg) analgesia as previously described (Bank et al., 2015) to monitor T b on line. T b was recorded every 3 min.
Additionally, three hamsters were culled in a normothermic state (T b 36.1 • C ± 0.7 • C) at the same ZT as control group (Figure 1, group 5). The brain was dissected from each hamster, frozen on dry ice and stored at −80 • C.

Isolation of Total RNA
Hypothalamic blocks were cut from frozen tissues between Bregma −0.20 and −2.70 mm, laterally at the hypothalamic sulci and dorsally 3-4 mm from the ventral surface. Tissue samples were homogenized in 500 µl TriFast using a micropestle. Total RNA was obtained using peqGOLD Trifast TM (Peqlab, Erlangen, Germany) according to the manufacturer's instructions. Total RNA was purified with the Crystal RNA MiniKit (Biolabproducts, Bebensee, Germany) including an on-column digestion with RNase-free DNase (Qiagen, Hilden, Germany). RNA integrity was proven by gel electrophoresis, total RNA was quantified spectrometrically and RNA purity was assessed by the 260/280 nm ratio on a NanoDrop 1000 spectrophotometer.

Illumina Sequencing
In total, 2 µg total RNA per sample were used for transcriptome analysis. Library preparation and Illumina sequencing were performed by GENterprise Genomics (Mainz, Germany). For library preparation the TruSeq RNA Library Preparation Kit (Illumina, San Diego, CA) was used. All RNA samples had a RIN ≥ 6.9. The samples were sequenced by Illumina NextSeq 500 with a calculated output of 50 million paired-end reads (2 × 150 bp) per sample. The raw Illumina data are available at the NCBI SRA database under the accession numbers biosample: SAMN062002211 to SAMN06200226 (Bioproject PRJNA360070). Since currently no annotated Phodopus sungorus genome is available, the reads were mapped against the genome of the Chinese hamster (Cricetulus griseus), that showed best compliance, using the CLC-Genomics Workbench 7.5.1 (Qiagen, Hilden, Germany). Only reads with intact pairs mapping with an 85% read identity and 85% read length were used for RPKM (reads per kilobase per million mapped reads) calculation (Mortazavi et al., 2008). Supplementary Table 1 shows the total number of reads and the number of reads mapped in pairs for each sample. Statistically significant expression changes between normothermic and torpid hamsters were calculated by an empirical analysis of digital gene expression (DGE) statistics, performing an "Exact Test" (Robinson and Smyth, 2008). This tool is implemented in CLC-Genomics Workbench 7.5.1. To correct for multiple testing, a false discovery rate (FDR) correction of p-values was applied (see Supplementary Table 2).
Transcripts with an RPKM-value ≥ 0.1 in one of the samples were chosen for further analysis. Transcripts with a fold change ≥ 1.2 and an FDR-corrected p ≤ 0.05 were considered as differentially expressed. The identified differentially expressed transcripts were functionally classified using the PANTHER Classification System (Protein ANalysis THrough Evolutionary Relationships; www.pantherdb.org) version 10.0 (Mi et al., 2013). Differentially expressed transcripts were additionally analyzed using the PANTHER Overrepresentation Test (release 20160715) applying the PANTHER GO-slim terms as annotated, followed by Bonferroni correction for multiple testing. The PANTHER Overrepresentation Test was conducted for all 284 differentially expressed genes as well as the 181 up regulated genes and the 103 down regulated genes respectively. Mus musculus was selected as reference organism for the GO annotation and for the statistical calculation of overrepresented GO-terms. Overrepresented terms with a Bonferroni corrected p ≤ 0.05 were considered as significant.

Experiment 2: Relative Quantification of Hypothalamic Gene Expression in Different Torpor Stages
Sampling To validate our Illumina results, we selected eight genes for further investigations by qPCR analysis. A group of seven genes immediately attracted attention for their potential role in structural changes (five collagens, myosin and dynein). Additionally, the von Willebrand factor (vwf ) was chosen for its role in blood clotting. To provide more detailed information about gene expression changes over a circadian cycle in torpid and normothermic state, 40 hamsters were used for gene expression analysis by real-time PCR (qPCR). A total of 20 animals were culled by carbon dioxide on a day with torpor expression at ZT1 (entrance into torpor: T b 30.8 • C ± 0.5 • C, n = 5), ZT4 (mid torpor: T b 22.5 • C ± 1.5 • C, n = 5), ZT7 (arousal: T b 30.4 • C ± 0.4 • C, n = 5) and ZT16 (active phase after torpor bout: T b 35.7 • C ± 0.6 • C, n = 5) (Figure 1, groups 1-4). Five normothermic animals were sampled for each time point as respective control group (ZT1: Tb 35.7 • C ± 0.5 • C; ZT4: (Figure 1, groups 5-8). Brains were dissected, frozen on dry ice and stored at −80 • C for qPCR analysis.

Isolation of RNA and cDNA Synthesis
Hypothalami were dissected from frozen brains and isolation of total RNA was performed as described for Experiment 1. For qPCR templates and generation of standard plasmids, cDNA was synthesized from every total RNA sample using RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA) and oligo-(dT)18 oligonucleotide primers following manufacturer's instructions. Reverse transcription was conducted using 1 µg total RNA per sample.

Real-Time qPCR and Analysis of Expression Data
qPCR was performed using Power SYBR R Green PCR Master Mix (Applied Biosystems, Darmstadt, Germany) on an ABI Prism 7300 Real Time PCR System (Applied Biosystems). Due to the large number of samples, qPCRs were performed on two 96-well plates for each target gene. For comparability, the normothermic ZT16 group was applied to all plates as interplate calibrator. Hprt was selected as reference gene, based on the stability of expression values across all samples. All samples were run in triplicates (for 5 biological replicates per group), using 1 µl cDNA as template in a reaction volume of 20 µl, and a series of six 10-fold dilutions of specific standard plasmids were used to generate the standard curve to calculate qPCR efficiencies. Additionally, a no-template control was included on each plate in duplicates for each target gene. Quantification was performed with the following cycling parameters for 40 cycles: 50 • C 2 min; 95 • C 10 min; 95 • C 15 s; 60 • C 15 s; 72 • C 30 s. Amplification specificity was controlled by dissociation curve analysis referring to the qPCR run.
First evaluation of qPCR results was carried out using the 7300 System Software v1.4.0 (ABI Prism, Applied Biosystems) and subsequently exported to Microsoft Excel 2010 to identify differences in expression levels using the CT method. All statistical testings and figures were done with SigmaPlot 12.5 (Systat Software Inc.). All results were statistically analyzed by two-way ANOVA with time of day (Zeitgebertime) and metabolic state (torpid/normothermic) as factors, followed by Tukey's test for pairwise comparison of relative expression levels between torpid and normothermic hamsters and within the torpid and normothermic groups.
The PANTHER overrepresentation test showed significant enrichments of the GO-slim terms only for the up regulated group of genes, comprising "transmembrane transporter activity" (9.64-fold, p = 0.034) in the domain molecular function and "extracellular matrix" (6.02-fold, p = 0.000397) in the domain cellular component.

Analysis of Most Affected Genes during Torpor Entrance
To determine the most affected genes during torpor entrance, we ranked the identified genes into the 20 most up-and 20 most down regulated genes, based on their fold changes.

Relative Gene Expression Patterns over the Circadian Cycle in Torpid and Normothermic Hamsters
To determine, whether differential candidate gene expression is restricted to torpor entrance and to assess circadian regulation, we investigated relative mRNA expression at ZT1, ZT4, ZT7, and ZT16 in animals undergoing torpor and animals remaining normothermic. Differences within each investigated time point are shown relative to normothermic control group at same ZT respectively (Figures 4A,C,E,G). Circadian variations for normothermic animals are shown relative to the normothermic ZT1 group. Circadian variations for torpid animals are presented relative to torpor ZT1 group (Figures 4B,D,F,H).
There was no effect of time of day on col17a1 mRNA levels for normothermic animals, but there was an effect of time of day for torpid animals (two-way ANOVA: p < 0.001). Col17a1 mRNA expression was reduced in the post-torpor group (ZT16) as compared to torpor entrance (ZT1, Tukey's test: p = 0.004) and to arousal (ZT7, Tukey's test: p = 0.002) (Figure 4B).
FIGURE 3 | Comparison of NGS-and qPCR data during torpor entrance. Expression changes were calculated by comparison of torpid animals at ZT1 to normothermic animals at ZT1 in both experiments. Gray bars represent torpid animals at ZT1(±SEM) analyzed by NGS and black bars represent torpid animals at ZT1(±SEM) analyzed by qPCR. Significant differences to their relative control groups are marked with *p < 0.05, **p < 0.01, ***p < 0.001.

DISCUSSION
Our data show 284 differentially expressed genes out of 27,830 identified genes in the hypothalamus of P. sungorus during entrance into the torpid state, implying that just a small set of genes is affected by the metabolic depression initiating torpor entrance. These results are in accordance with previous studies showing that transcript levels of most genes are unaffected during torpor (Storey and Storey, 2004). In accordance with the fact that daily torpor is a state of extreme metabolic adjustment, the majority of differentially regulated genes was found in cellular and metabolic processes for both, up and down regulated genes.
The majority of the top 20 down regulated genes were transcription factors, which could be responsible for a delay or suppression of mRNA transcription during the torpid state. It has been shown before, that transcriptional initiation as well as elongation rates are reduced during hibernation in goldenmantled ground squirrels (van Breukelen and Martin, 2002). Also in P. sungorus metabolic depression is associated with reduced transcriptional initiation (Berriel Diaz et al., 2004). This may contribute to the generally suppressed protein synthesis during torpor that has been demonstrated in various tissues from different species (Gulevsky et al., 1992;Frerichs et al., 1998;Hittel and Storey, 2002).
Within the top 20 up regulated group our data show a remarkable number of genes coding for structure proteins. Except for the up regulation in collagen genes we were able to verify these results by qPCR for dnah2, myo15a and the procoagulation factor vwf.
Collagens are extracellular matrix structural components, which are involved in neuronal development of the brain.
Collagens play a role in axonal guidance, synaptogenesis and establishment of brain architecture (Chernousov et al., 2006;Fox, 2008;Hubert et al., 2009). A study of Schwartz et al. (2013) identified an up regulation of several collagen genes in the cerebral cortex, but not in the hypothalamus, of thirteen-lined ground squirrels during deep hibernation and interbout arousals, indicating synaptic plasticity during hibernation (Schwartz et al., 2013). Although we obtained a significant up regulation in mRNA expression of five collagen genes during torpor entrance by NGS and a significant enrichment of extracellular matrix components in up regulated gens, we were not able to verify these results by qPCR. There was only a trend of increased col17a1 during torpor entrance as well as slightly lower mRNA levels during all other torpor stages and no diurnal changes could be detected in normothermic animals. Investigation of all other collagens identified in the 20 up regulated group showed a similar picture with a slight up regulation at torpor entrance and trend to lower mRNA levels during the other torpor stages that did not reach significance (data not shown). There was a high variability in the mRNA expression levels of qPCR samples, especially at torpor entrance, which might have caused the non-significant result. Different groups of animals were used for NGS and qPCR study and data might reflect inter-individual differences. A larger sample size might help to resolve expression patterns in collagen genes more precisely. Hence, whether collagens are involved in synaptic remodeling and plasticity during torpid states remains to be revealed.
Elevated expression of dnah2-and myo15a mRNA during torpor entrance could be identified by both, NGS and qPCR approach. Myosin and dynein are structural components of cytoskeleton and represent two out of three superfamilies of molecular motor proteins in neurons. They are able to transport biomolecules, such as vesicles, protein complexes and mRNAs in axons, dendrites and pre-and post-synaptic regions. Intracellular transport is necessary for neuronal morphogenesis, function and survival (Hirokawa et al., 1998(Hirokawa et al., , 2010Vale, 2003). During deep hibernation, elevated mRNA levels of three different myosin types and one dynein have been detected in the cerebral cortex of S. tridecemlineatus, indicating dynamic structural changes (Schwartz et al., 2013).
In our study, hamsters showed elevated dnah2 and myo15a expression only during torpor entrance (ZT1), whereas mRNA expression was reduced at all other investigated torpor states (mid torpor, arousal, post-torpor) compared to normothermic animals. The higher expression of dnah2 and myo15a during torpor entrance could be important to ensure maintenance of synaptic transmission and neuron survival during torpor by an elevated transport of biomolecules. It might also be possible that higher mRNA amounts are produced at the beginning and stored during the torpid state to provide transcripts for a fast utilization of these molecular motors during arousal. However, we think this possibility is unlikely because mRNA levels are already declining during mid-torpor (ZT4).
In normothermic animals dnah2 as well as myo15a show a diurnal regulation in its mRNA expression with a peak at ZT7 in normothermic animals. This might suggest a higher demand of these motor proteins during the hamster's naturally active phase.
FIGURE 4 | Circadian regulation of col17a1, dnah2, myo15a, and vwf in torpid and normothermic Djungarian hamsters. Bar graphs on the left side show differences in mRNA expression of col17a1 (A), dnah2 (C), myo15a (E), and vwf (G) in torpid animals (gray bars, ±SEM) relative to normothermic control group at same ZT (black bars, ±SEM). Significant differences are marked with *p < 0.05, **p < 0.01, and ***p < 0.001. Line graphs on the right side show relative differences in mRNA expression of col17a1 (B), dnah2 (D), myo15a (F), and vwf (H) over the course of a day within normothermic animals relative to normothermic ZT1 group marked with upper case (black circles, ±SEM) and within torpid animals relative to torpid ZT1 group marked with lower case (gray circles, ±SEM). Data points with different characters are significantly different (p < 0.05).
Taken together changes in structural protein shows evidence for plasticity in the hypothalamus of torpid hamsters and thereby confirm studies in deep hibernation that have proposed plastic changes in the brain before.
In addition to structure gene expression changes, we chose to investigate vwf in more detail, because of its function in blood clotting. In torpid animals the reduced heart rate, ventilation and T b results in a decreased blood flow that increases relatively fast to its euthermic flow rate during arousal. In contrast to all other mammalian species, torpor expressing mammals are able to survive these periods of low blood flow and consequent reperfusion without apparent formations of deep vein thrombi, stroke or pulmonary embolism (Lyman and O'Brien, 1961;Frerichs et al., 1994).
vWF is a major factor involved in platelet adhesion and thrombus formation (Denis and Wagner, 2007). Higher vWF levels increase the risk for thrombosis and embolism whereas deficiency in vWF activity leads to the human bleeding disorder von Willebrand's disease (Sadler, 1998(Sadler, , 2005. Moreover, Zhao et al. (2009) identified vWF as an important protein regulating the occurrence of cerebral ischemia and showed that a lack of vWF is able to reduce infarct volume (Zhao et al., 2009). Based on this knowledge, a reduced level of vWF would be expected during the torpid state to prevent blood clotting during periods of low blood flow. Indeed, in plasma samples of hibernating thirteen-lined ground squirrels vWF collagen binding is 10fold decreased and in lung tissues vwf mRNA expression is 3-fold down regulated during torpor (Cooper et al., 2016). Unexpectedly, our NGS and qPCR data show an elevated level of vwf mRNA during torpor entrance in the hypothalamus. The elevated level of vwf mRNA might either not directly translate into protein variation or alternatively translate into protein without damaging effects, namely inactive vWF. vWF is a large multimeric glycoprotein which can be cleaved in smaller multimers by ADAMTS13, a zinc-containing metalloprotease enzyme. These smaller multimers of vWF have a strongly decreased activity resulting in a reduced platelet adhesion and aggregation (Chauhan et al., 2006;Zhao et al., 2009). In this case, no damage of brain structures would be expected even when higher vwf levels are present. Moreover, apart from the up regulation during torpor entrance, vwf expression was lower in torpid animals at all other investigated states, supporting the hypothesis of low vwf levels facilitating blood flow during the torpid state. Diurnal changes of vwf could be detected in either group. Normothermic animals displayed highest vwf level at ZT7, suggesting a higher demand of vwf at the beginning of the active time. In torpid animals vwf level is lowest at mid torpor (ZT4) and post-torpor (ZT16). Taken together, our data provide evidence for readjustment of blood clotting during different torpor stages as well as times of day.
In general, the diurnal mRNA expression of all investigated genes of this study is less pronounced in torpid animals, which is likely to be caused by the suppression of transcription and translation during torpor. The transcriptional depression during torpor has been shown to result from both, down regulated transcriptional initiation and suppressed elongation (van Breukelen and Martin, 2002;Berriel Diaz et al., 2004). Low T b during torpor affects biochemical process, leading to a decline in gene expression caused by the temperature sensitivity of transcriptional elongation (van Breukelen and Martin, 2002;Berriel Diaz et al., 2004).
The NGS technology allows a whole transcriptome survey of gene expression changes and our analysis provide an overview of gene expression changes during torpor initiation in P. sungorus for the first time.
Although we could not determine signaling pathways regulating torpor initiation with this approach, we identified molecular adaptions in the hypothalamus of P. sungorus initiated during the early state of torpor. Our data provide evidence for synaptic remodeling and plasticity, an elevated transport of biomolecules and readjustment of coagulation. Comparable gene expression changes have already been found in deep hibernators. This would support the hypothesis that daily torpor and hibernation are similar physiological states only differing in amplitude and duration. Interestingly, the molecular changes already occur within the short time span of daily torpor. These adaptations may, just like in deep hibernation, help the brain cells to better survive or reduce cell damages during the extreme physiological conditions in the torpid state. In the future, precise anatomical investigation of identified genes is necessary to eventually gain insights into their functions.

AUTHOR CONTRIBUTIONS
CC, AF, and AH designed experiments, CC and JK performed experiments. CC, JK, AF, and AH analyzed and interpreted the data. CC and AH drafted the manuscript which was critically revised by JK and AF.

FUNDING
This work was funded by the German Research Foundation (DFG, Emmy-Noether HE6383 to AH).