the expression of Per3 and Rorα in any ZT, while higher gene expression were observed for Clock at ZT3 compared with ZT11 (P = 0.02), for Bmal1 at ZT17 compared with ZT11 (P = 0.01), for Cry1, and Cry2 at ZT17 compared with ZT3 (P = 0.03 and P = 0.02, respectively), and for Per1 at ZT11 compared with ZT3 (P = 0.02). Regarding Per2, the highest expression was observed at ZT11 and the lowest at ZT17 in all groups (P < 0.01). The expression of Clock, Per1, and Per3 genes did not show differences at ZT3, ZT11, and ZT17 in CG; however, in RF-d and RF-n groups, these genes presented higher expression in the morning and in the afternoon (P < 0.01). Compared with CG and RF-n groups, in which food access occurred without dissociation with the usual pattern of rat nocturnal activity, higher expression of Bmal1 was observed at ZT11 (P = 0.05) in the RF-d group, while Cry1 (P < 0.01), Cry2 (P = 0.01), and Rorα at ZT3 and ZT11 (P < 0.01).

The expression of clock genes in the PVN of CG, RF-n, RF-d, and DF groups at ZT3, ZT11, and ZT17 are presented in **Figure 4**. Regarding differences among morning, afternoon, and night values in the same feeding schedule, we observed no differences in the expression of Clock, Bmal1, and Cry1 genes at ZT3, ZT11, and ZT17 in any studied group. In the CG, higher Per1 expression was observed in the afternoon (ZT11) and at night (ZT17) (P < 0.01), while higher expression of Per2, Per3, and Rorα (P < 0.01) were observed at night (ZT17). No different Rev-erbα expression was observed among the three different ZT points. In the RF-n group, the pattern of expression of Per2 (P < 0.01), Per3 (P = 0.04), and Rorα (P < 0.01) genes was similar to CG, with higher expression at night (ZT17), while Per1 and Cry2 diurnal variation observed in the CG was abolished in this group. In the RF-d group, there was no difference in the expression of the studied genes among ZT3, ZT11, and ZT17, with exception of Rev-erbα, which had higher expression at night (ZT17) (P = 0.03). DF group showed higher expression of Rorα and Rev-erbα at ZT3 compared with ZT17 (P = 0.02).

The expression of clock genes in the ARC of CG, RF-n, RFd, and DF groups at ZT3, ZT11, and ZT17 are presented in **Figure 5**. No difference in Bmal1 expression was observed among ZT3, ZT11, and ZT17 in any group. In the CG, Cry1 expression was higher (P = 0.02) in the morning (ZT3), while Clock, Per1, Per2, Per3, Cry2, Rev-erbα, and Rorα (P < 0.01) showed higher expression at night (ZT17).  RF-n also had higher Per1 expression at ZT17 (P < 0.01), while RF-d exhibited higher Per1 expression at morning (ZT3) and at night (ZT17) (P < 0.01). RF-d group maintained the pattern of Rev-erbα expression observed in CG group, while RF-n lost this pattern. DF group presented higher expression of Cry1 at ZT17 than at ZT3 and ZT11 (P < 0.01).

#### The Feeding Restriction Schedules and the Modulation of Temporal Relationship of Clock Gene Expression with Corticosterone Secretion

In the SCN of control animals, the higher Per2 expression was concomitant with the higher corticosterone secretion while no relationship was observed regarding Clock, Bmal1, Per1, Per3, Cry1, Cry2, Rev-erbα, and Rorα genes. In the RF-n group, besides Per2, the expression of Clock, Per1, and Per3 genes also exhibited higher expression coinciding with higher corticosterone secretion. The same pattern was observed for Per1 and Per3 in RF-d group, in which there was also a delay of the higher expression of Clock, Bmal1, Per2, Cry1, Cry2, and Rorα compared with the higher corticosterone secretion. On the other hand, in DF group the most interesting finding was the anticipatory higher expression of Clock and Bmal1 regarding the corticosterone secretion.

In the PVN of control animals, the higher expression of Per1, Per2, Per3, Cry2, and Rorα were concomitant with the higher corticosterone secretion, while no relationship was observed regarding Clock, Bmal1, Cry1, and Rev-erbα genes. In the RFn, the higher expression of Per2, Per3, and Rorα was delayed concerning the corticosterone secretion, while there was no relationship regarding the expression of Clock, Bmal1, Per1, Cry1, Cry2, and Rev-erbα genes. On the other hand, in the RF-d and DF groups, the majority of the studied genes lost

their relationship with corticosterone secretion compared with the pattern observed in control and RF-n groups. Of note, in DF group, the higher expression of Rorα and Rev-erbα was anticipatory concerning corticosterone secretion.

In the ARC of control animals, while the higher expression of Clock, Per1, Per2, Per3, Cry2, Rev-erbα, and Rorα was coincident, the Cry1 higher expression was anticipatory in relation to corticosterone secretion. No relationship was observed regarding Bmal1 gene. On the other hand, in the RF-n and RF-d groups, there was a delay of the higher expression of Clock, Per1, Per3, Cry2, Rev-erbα, and Rorα compared with corticosterone secretion. No relationship was observed concerning Bmal1, Cry1, and Per2 in RF-n group and Rev-erbα in RF-d group, respectively. Of note, in DF group, most of the genes lost the pattern of expression, also missing any relationship with corticosterone secretion.

#### DISCUSSION

In the present study, rats submitted to three different schedules of food restriction changed the diurnal pattern of Per1, Per2, and Per3 expression in the SCN and in the extra-SCN oscillators— PVN and ARC nuclei—suggesting that non-photic cues can modulate the expression of clock genes, but not in a uniform way, throughout the hypothalamic nuclei. Furthermore, the corticosterone anticipatory peak was verified irrespective to the food timing access, except in the DF group. Altogether our data confirm that restricted feeding schedules can be an entrainment cue for glucocorticoid circadian variation.

The initial body weight of rats was similar in all groups, but at the end of the experiment, RF-d and RF-n groups showed lower weight as well as lower daily food intake than control and DF groups, similarly to studies using food restriction but with different duration of food access (Honma et al., 1983; Stephan and Becker, 1989; Stephan, 2002). In our study, the time of feeding in RF-d and RF-n groups was restricted by 2 h of food access leading to about 50% calorie intake restriction when compared with controls, suggesting that body weight could be determined primarily by the amount of energy intake and not by the time of feeding. Indeed, control and DF groups showed similar body weight. Arble et al. (2009) showed that mice fed during the 12-h light phase gained even more weight than mice fed on an equivalent amount of calories during the 12-h dark

\*P < 0.05, \*\*P < 0.01.

phase (Arble et al., 2009). Of note, in that study, differently from our study, mice were fed with high-fat diet during 6 weeks. Altogether, these data point out that a longer extent of daytime feeding would be necessary to achieve higher body weight.

Change in the daily eating pattern is considered a powerful Zeitgeber for the diurnal rhythm of glucocorticoids (Krieger, 1974; Leal and Moreira, 1996; Girotti et al., 2009). Higher corticosterone concentrations were observed in the control and RF-n groups at ZT11 than at ZT3, coincident with the onset of the nocturnal activity period of rats (Le Minh et al., 2001). On the other hand, RF-d group, which was allowed to eat only during 2 h in the morning, showed higher concentrations of corticosterone at ZT3 than at ZT11, inverting the rat corticosterone circadian pattern, as previously described (Krieger, 1974; Moreira and Krieger, 1982; Leal et al., 1995; Leal and Moreira, 1996). Of note, in the RF-n group, in which food access was restricted but occurred without dissociation with the pattern of rat nocturnal activity, we observed an entrainment of corticosterone peak coincident with food availability. Furthermore, the ZT3 and ZT11 time points allowed us to observe the well known anticipatory corticosterone peak before food availability during the light phase of the light/dark cycle (Oliveira et al., 1993; Le Minh et al., 2001). This finding was not observed in the DF group, which presented higher corticosterone levels at ZT11 as control and RF-n groups. Our data are in agreement with previous reports showing that anticipatory corticosterone peak did not happen when food availability occurs during longer intervals of the day (Honma et al., 1983; Belda et al., 2005). Probably, in this situation the photic cue continues as the main Zeitgeber.

In the SCN, all groups of restricted feeding schedules exhibited different patterns of the expression of Clock and Per1 genes when compared to the control group, while the Per3 gene expression was only modified in the RF-d and RF-n groups, characterized by an intense food deprivation (food availability restricted by 2 h). The expression of Bmal1, Cry1, and Cry2 was only modified in the RF-d and DF groups, which food access occurred in dissociation with the pattern of rat nocturnal activity, presenting higher expression in the morning, and early evening. In contrast, a higher expression of Per2 at ZT11 than at ZT17 was observed in all groups, thus the expression of this gene was

unchanged by variable food schedule. Our data demonstrated, in agreement with Caldelas et al. (2005) and Minana-Solis et al. (2009) but in contrast to others (Damiola et al., 2000; Hara et al., 2001; Wakamatsu et al., 2001; Verwey et al., 2007) that SCN is also affected by food entrainment. This controversial data can be explained mainly because of the period in which rats were food entrained in the different protocols, as well as by the duration of daily food access, which varied from 2 h, similar to our study, as long as 12 h as established by Damiola et al. (2000). Thus, time of the day as well as the duration of daily food access may modulate the SCN response, as previously suggested by Minana-Solis et al. (2009). Therefore, feeding cycle can be considered the dominant Zeitgeber for the clock gene expression not only in the peripheral but also in the central pacemakers.

We were not able to detect daily variations in the expression of clock genes in the SCN of the control group, except for Per2. Our protocol explored three time points only, which could be a possible limitation to properly study variations in circadian rhythm but it is still informative to understand the role of restricted feeding on the expression of clock genes in hypothalamic nuclei. However, in the PVN and ARC, higher expressions of Per2, Per3, and Cry2 were observed synchronically at ZT17 in the dark period, whereas the other clock genes phased differently. In addition, different food schedules changed the expression of some clock genes in the SCN, PVN, and ARC, compared with CG. Thus, our data demonstrated that under ad libitum food access or under food restriction there are different phases in the clock gene expression between the SCN and extra-SCN nuclei as PVN and ARC hypothalamic nuclei.

Data on ARC and in other energy balance and feeding behavior related areas have shown that dietary restriction induces or increases the amplitude of the daily rhythms of Per1 and Per2 genes and their proteins (Mieda et al., 2006; Verwey et al., 2007; Angeles-Castellanos et al., 2008; Guilding et al., 2009). In addition, c-fos, Per1, Per2, and Bmal1 rhythms were in antiphase with the respective SCN genes in the PVN in rats under food restriction, suggesting that extra-SCN neuronal clocks can be, but not always, expressed in antiphase from the SCN (Damiola et al., 2000). Another study also showed that the daily food restriction did not lead to a uniform timing effect in clock gene expression in the different hypothalamic nuclei (Minana-Solis et al., 2009). Our data expand these findings since we studied a larger number of clock genes and three different food restriction schedules. Indeed, RF-n, RF-d, and DF groups also showed differential effects on the expression of clock genes, inducing phase alteration in some cases or keeping them unaffected in others. Of note, daily restricted feeding can produce food-entrainable oscillations in the cerebral cortex, hippocampus, striatum, pyriform cortex, and PVN even in SCN-lesioned animals (Wakamatsu et al., 2001). Thus, clock genes in extra-SCN oscillators do not respond to feeding schedules in a uniform manner, suggesting a tissuespecific regulation (Feillet et al., 2008; Minana-Solis et al., 2009).

Our data regarding the modulation of feeding restriction schedules on temporal relationship between clock gene expression and corticosterone secretion show that the most important changes were observed in RF-d and DF groups, in which conditions the eating pattern was not coincident with the onset of the nocturnal activity period, suggesting that a disruption in feeding cues may cause dissociation between gene expression and corticosterone secretion. This misalignment can explain night eating syndrome observed in humans, which is characterized by increased food consumption late in the day, in the evening, and in the night and has been linked to obesity and increased risk of metabolic syndrome (Nader et al., 2010).

In conclusion, our data show that feeding restriction schedules induce changes in the body weight, in the corticosterone circadian variation, and in the expression of the clock genes in hypothalamic nuclei, suggesting the ability of animals to predict, anticipate, and respond to food availability. The pattern of clock gene expression in response to RF-n, RF-d, and DF schedules was not affected uniformly, suggesting that the master oscillator in SCN as well as PVN and ARC, areas involved in energy homeostasis and homeostatic aspect of feeding, are regulated in a tissue-specific manner in response to feeding restriction, a non-photic cue, in order to entrain clock gene system.

#### REFERENCES

### AUTHOR CONTRIBUTIONS

LD: Conception and design of research, performed experiments, analyzed data, interpreted results of experiments, prepared figures, drafted manuscript, edited, and revised manuscript, approved final version of manuscript. SR: Performed experiments, interpreted results of experiments, prepared figures, revised manuscript and approved final version of manuscript. AB: Performed experiments, analyzed data and interpreted results of experiments. FC: Performed experiments, analyzed data and interpreted results of experiments. CM: Performed experiments, analyzed data and interpreted results of experiments, revised manuscript and approved final version of manuscript. EU: Conception and design of research, performed experiments and interpreted results of experiments. JA: Conception and design of research, revised manuscript, approved final version of manuscript. LE: Conception and design of research, revised manuscript, approved final version of manuscript. PE: Conception and design of research, revised manuscript, approved final version of manuscript. AM: Conception and design of research, analyzed data, interpreted results of experiments, edited and revised manuscript, approved final version of manuscript. MC: Conception and design of research, analyzed data, interpreted results of experiments, prepared figures, drafted manuscript, edited and revised manuscript, approved final version of manuscript. I certify that all authors had a substantial contribution to the manuscript.

#### FUNDING

National Counsel of Technological and Scientific Development— CNPq, Coordination for the Improvement of Higher Level Personnel—CAPES, Sao Paulo Research Foundation—FAPESP (Process numbers: 2007/58365-3, 2013/09799-1).

#### ACKNOWLEDGMENTS

We thank Rogerio Zuliani, Jose Roberto Silva, Maria Valci dos Santos, and Rubens Fernando de Melo for their excellent technical assistance.

hormones and the autonomic nervous system. J. Endocrinol. 177, 17–26. doi: 10.1677/joe.0.1770017

with insulin in daily energy flow. Front. Neuroendocrinol. 14, 303–347. doi: 10.1006/frne.1993.1010

circadian pacemaker. Proc. Natl. Acad. Sci. U.S.A. 103, 12150–12155. doi: 10.1073/pnas.0604189103

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 De Araujo, Roa, Bueno, Coeli-Lacchini, Martins, Uchoa, Antunes-Rodrigues, Elias, Elias, Moreira and De Castro. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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

Ceyda Cubuk, Julia Kemmling, Andrej Fabrizius and Annika Herwig\*

Zoologisches Institut, Universität Hamburg, Hamburg, Germany

#### Edited by:

Etienne Challet, CNRS UPR3212, University of Strasbourg, France

#### Reviewed by:

Fabien Pifferi, UMR 7179 Centre National de la Recherche Scientifique, Muséum National d'Histoire Naturelle, France Dietmar Weinert, Martin Luther University of Halle-Wittenberg, Germany

> \*Correspondence: Annika Herwig<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Neuroendocrine Science, a section of the journal Frontiers in Neuroscience

Received: 28 October 2016 Accepted: 27 February 2017 Published: 13 March 2017

#### Citation:

Cubuk C, Kemmling J, Fabrizius A and Herwig A (2017) Transcriptome Analysis of Hypothalamic Gene Expression during Daily Torpor in Djungarian Hamsters (Phodopus sungorus). Front. Neurosci. 11:122. doi: 10.3389/fnins.2017.00122 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.

Keywords: metabolic depression, seasonal adaptations, circadian, hypothalamus, Illumina, RNA-Seq

**Abbreviations:** col17a1, collagen, type XVII, alpha 1; dnah2, dynein, axonemal, heavy chain 2; myo15a, myosin XVA; NGS, next generation sequencing; qPCR, quantitative real-time PCR; Ta, ambient temperature; Tb, body temperature; vwf, von Willebrand factor; ZT, Zeitgebertime.

## 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 (Ta) 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<sup>a</sup> 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 (Tb) 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 (Kirsch et al., 1991). 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 (Ouarour et al., 1991; 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., 2011, 2013; Schwartz et al., 2013, 2015; Grabek 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.

### MATERIALS AND METHODS

#### 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 Ta. 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<sup>a</sup> 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<sup>b</sup> on line. T<sup>b</sup> was recorded every 3 min.

## Experiment 1: Transcriptomic Analysis of Hypothalamic Gene Expression at Torpor Entrance

#### Sampling

Between week 13 and week 16, three animals were euthanized by carbon dioxide during entrance into torpor (T<sup>b</sup> 30.4◦C ± 0.6◦C) at Zeitgebertime 1 (ZT1; ZT0 = lights on) (**Figure 1**, group 1).

Additionally, three hamsters were culled in a normothermic state (T<sup>b</sup> 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 TrifastTM (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<sup>b</sup> 30.8◦C ± 0.5◦C, n = 5), ZT4 (mid torpor: T<sup>b</sup> 22.5◦C ± 1.5◦C, n = 5), ZT7 (arousal: T<sup>b</sup> 30.4◦C ± 0.4◦C, n = 5) and ZT16 (active phase after torpor bout: T<sup>b</sup> 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: T<sup>b</sup> 35.7◦C ± 0.4◦C; ZT7: T<sup>b</sup> 35.6◦C ± 0.4◦C;

ZT16: T<sup>b</sup> 36.2◦C ± 1.3◦C) (**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.

#### Cloning and Sequencing

For standard plasmids, coding sequence fragments (100–200 bp long) of collagen alpha-1(XXIV) chain-like (LOC103164493), collagen, type XX, alpha 1 (col20a1), collagen, type XVII, alpha 1 (col17a1), collagen, type XVIII, alpha 1 (col18a1), collagen, type V, alpha 3 (col5a3), dynein, axonemal, heavy chain 2 (dnah2), myosin XVA (myo15a), von Willebrand factor (vwf) and hypoxanthine phosphoribosyltransferase (hprt) were amplified by gene specific primers (**Table 1**). All primers were designed on the P. sungorus specific sequences generated by Illumina sequencing using the onlinetool OligoAnalyzer 3.1. The primers were designed with a melting temperature at 60◦C ± 1.1◦C. After cloning of the amplicons using the pGEM <sup>R</sup> -T Easy Vector System (Promega, Madison, USA) following the manufacturer's instructions, the cDNA fragments were Sanger-sequenced by GATC Biotech (Konstanz, Germany).

#### TABLE 1 | P. sungorus specific primer sequences.

#### Real-Time qPCR and Analysis of Expression Data

qPCR was performed using Power SYBR <sup>R</sup> 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 11CT 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.

## RESULTS

### Transcriptomic Expression Analysis in the Hypothalamus of P. sungorus during Torpor Entrance

We identified a total number of 27,830 transcripts with 284 transcripts being differentially expressed in hamsters during torpor entrance as compared to ZT matched normothermic hamsters (**Table 2**). A total of 181 transcripts were significantly upregulated whereas 103 transcripts were significantly downregulated. All transcripts identified had an RPKM-value ≥ 0.1 and a FDR p ≤ 0.05.

### Functional Classification of Differentially Expressed Genes

Differentially expressed transcripts during torpor entrance were classified according to gene ontology categories. It has to be noted, that some genes were included in more than one category.

Of the 181 up regulated genes 154 could be classified and assigned to 389 biological processes (**Figure 2A**). The majority of up regulated genes are involved in metabolic (18%) and cellular processes (18%), followed by localization (12%), biological regulation (9%), developmental processes (9%), cellular component organization or biogenesis (7%), multicellular organismal processes (7%), biological adhesion (6%), immune system processes (5%), response to stimulus (5%), apoptotic processes (2%), and reproduction (2%).

TABLE 2 | Overview of Illumina sequencing data and transcriptomic expression analysis.

A total number of 82 out of the 103 down regulated genes could be classified and assigned to 117 biological process hits (**Figure 2B**). The majority of down regulated genes was assigned to metabolic processes (33%), followed by cellular processes (20%), biological regulation (10%), developmental processes (9%), multicellular organismal processes (5%), localization (4%), immune system processes (4%), cellular component organization or biogenesis (3%), reproduction (3%), response to stimulus (3%), apoptotic processes (3%), and biological adhesion (3%).

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.

Most up regulated genes (**Table 3**) showed fold changes in a range of 1.55–2.66. Within this group we found 8 genes coding for structure proteins (LOC103164493, col20a1, myo15a, col17a1, micalcl, dnah2, col18a1, col5a3), 4 involved in transporter function (abca6, atp2a1, kcnh3, atp1a4), 2 with signaling function (OR2K2, LOC100773864) and one gene each involved in stress defense (klk8), coagulation (vwf) and cell death (steap3). Three genes have so far unknown function (LOC103160902\_1, LOC100766933, catip).

Fold changes of the most down regulated genes (**Table 4**) ranged between −1.5 and −4.0. This group contained 7 genes coding for transcription factors (stk31, LOC100768314, LOC102638674, LOC100756005, LOC102632383, LOC102642077, smim11), 4 with enzymatic activity (top2a, clk1, coq3, LOC100772408\_2), 2 with transporter

#### TABLE 3 | Most up regulated genes in the hypothalamus during torpor entrance.

functions (slc47a1, nipsnap3b), one gene each involved in cellular structure (cornifin-A), signaling (psmc3ip), rRNA maturation (rrp15) and ORF (swt1) and 3 genes with so far unknown function (LOC100754037, LOC103159055, LOC100753290).

#### Verification of Most up Regulated Gene Expression by qPCR

To verify the Illumina results, we calculated hypothalamic relative mRNA expression of col17a1, dnah2, myo15a, and vwf during torpor entrance by qPCR (**Figure 3**). Up regulation could be confirmed for dnah2 (qPCR: 1.6-fold, p = 0.016; Illumina: 1.7 fold, p < 0.001), myo15a (qPCR: 2.5-fold, p = 0.005; Illumina: 1.9-fold, p = 0.035) and vwf (qPCR: 1.6-fold, p < 0.001; Illumina: 1.6-fold; p = 0.046). Up regulation of col17a1 did not reach significance in the qPCR analysis (qPCR: 1.2-fold, p = 0.462; Illumina: 1.8-fold, p = 0.027). Also the other collagens identified by Illumina did not reach significance by qPCR (col5a3: 1.3-fold, p = 0.550; col18a1: 1.2-fold, p = 0.450; col20a1: 2.0-fold, p = 0.361; LOC103164493: 1.4-fold, p = 0.563) (data not shown).

### 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**).

There were no significant changes in mRNA expression during torpor entrance (ZT1), mid torpor (ZT4) or arousal (ZT7) relative to normothermic control groups at the same ZTs. Posttorpor (ZT16), mRNA expression was 0.57-fold down regulated (Tukey's test: p = 0.014) (**Figure 4A**).

There was an effect of time of day on dnah2 mRNA levels both, in normothermic and torpid animals (two-way ANOVA: p < 0.001). Normothermic animals showed lowest mRNA expression at ZT1 (Tukey's test: ZT1 vs. ZT4 p = 0.006, ZT1 vs. ZT7 p < 0.001, ZT1 vs. ZT16 p = 0.004) that increased at ZT4, peaked at ZT7 (Tukey's test: ZT7 vs. ZT4 p = 0.011) and decreased again at ZT16 (Tukey's test: ZT16 vs. ZT7 p = 0.015) (**Figure 4D**). Over the investigated torpor stages (ZT1, 4, 7) no significant changes were found, but post-torpor (ZT16) mRNA expression

TABLE 4 | Most down regulated genes in the hypothalamus during torpor entrance.

was 0.38-fold down regulated as compared to torpor entrance (ZT1, Tukey's test p < 0.001), mid torpor (ZT4, Tukey's test p = 0.003) and arousal (ZT7, Tukey's test p < 0.001) (**Figure 4D**).

Relative to their normothermic control groups, dnah2 expression was 1.64-fold up regulated during torpor entrance (ZT1, Tukey's test p = 0.016), 0.56-fold down regulated during mid-torpor (ZT4, Tukey's test p = 0.014), 0.40-fold down regulated during arousal (ZT7, Tukey's test p < 0.001) and 0.26 fold down regulated post-torpor (ZT16, Tukey's test p < 0.001) (**Figure 4C**).

There was an effect of time of day on myo15a mRNA levels both, in normothermic and torpid animals (two-way ANOVA: p < 0.001). Normothermic animals showed low mRNA expression at ZT1 that increased at ZT4 (Tukey's test p = 0.019) and ZT7 (Tukey's test p < 0.001), before decreasing again at ZT16 (Tukey's test: ZT16 vs. ZT7 p = 0.030). Torpid animals showed highest mRNA expression during torpor entrance (ZT1) differing significantly from mid torpor (ZT4, Tukey's test p = 0.010) and post-torpor (ZT16, Tukey's test p < 0.001). mRNA expression at mid torpor (ZT4) was also down regulated as compared to arousal (ZT7, Tukey's test p = 0.047) and mRNA expression during arousal (ZT7) was up regulated relative to post-torpor (ZT16, Tukey's test p < 0.001) (**Figure 4F**).

Myo15a expression was 2.51-fold up regulated during torpor entrance (ZT1, Tukey's test p = 0.005), 0.32-fold down regulated at mid torpor (ZT4, Tukey's test p = 0.002), 0.40-fold down regulated during arousal (ZT7, Tukey's test p = 0.006) and 0.29-fold down regulated in post-torpor group (ZT16, Tukey's test p < 0.001) (**Figure 4E**) as compared to the normothermic control groups.

There was an effect of time of day on vwf mRNA levels both, in normothermic and torpid animals (two-way ANOVA: p < 0.001). Normothermic animals showed low mRNA expression at ZT1 that increased at ZT4 (Tukey's test p = 0.038) and further at ZT7 (Tukey's test: ZT1 vs. ZT7 p < 0.001, ZT4 vs. ZT7 p = 0.017), before decreasing at ZT16 (Tukey's test: ZT7 vs. ZT16 p < 0.001). The mRNA expression in torpid animals was significantly up regulated during torpor entrance (ZT1) as compared to mid torpor (ZT4, Tukey's test p = 0.002) and post-torpor (ZT16, Tukey's test p < 0.001) and during arousal (ZT7) compared to post-torpor (ZT16, Tukey's test p = 0.002) (**Figure 4H**).

Vwf expression was 1.61-fold up regulated during torpor entrance (ZT1, Tukey's test p < 0.001), 0.58-fold down regulated during mid-torpor (ZT4, Tukey's test p = 0.004), 0.51-fold down regulated during arousal (ZT7, Tukey's test p < 0.001) and 0.60 fold down regulated in post-torpor group (ZT16, Tukey's test p = 0.003) (**Figure 4G**) as compared to the normothermic control groups.

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, 2010; Vale, 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.

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<sup>b</sup> 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, 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 10 fold 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.

#### REFERENCES

Bank, J. H., Kemmling, J., Rijntjes, E., Wirth, E. K., and Herwig, A. (2015). Thyroid hormone status affects expression of daily torpor and gene transcription in Djungarian hamsters (Phodopus sungorus). Horm. Behav. 75, 120–129. doi: 10.1016/j.yhbeh.2015.09.006

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<sup>b</sup> 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).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins. 2017.00122/full#supplementary-material

of ADAMTS13. J. Exp. Med. 203, 767–776. doi: 10.1084/jem.200 51732

and daily torpor in Djungarian hamsters (Phodopus sungorus). J. Endocrinol. 106, 367–375. doi: 10.1677/joe.0.1060367

Zhao, B. Q., Chauhan, A. K., Canault, M., Patten, I. S., Yang, J. J., Dockal, M., et al. (2009). von Willebrand factor-cleaving protease ADAMTS13 reduces ischemic brain injury in experimental stroke. Blood 114, 3329–3334. doi: 10.1182/blood-2009-03-213264

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Cubuk, Kemmling, Fabrizius and Herwig. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Role of Circadian Rhythms in Muscular and Osseous Physiology and Their Regulation by Nutrition and Exercise

#### Shinya Aoyama1, 2 and Shigenobu Shibata<sup>2</sup> \*

<sup>1</sup> Organization for University Research Initiatives, Waseda University, Tokyo, Japan, <sup>2</sup> Laboratory of Physiology and Pharmacology, School of Advanced Science and Engineering, Waseda University, Tokyo, Japan

The mammalian circadian clock regulates the day and night cycles of various physiological functions. The circadian clock system consists of a central clock in the suprachiasmatic nucleus (SCN) of the hypothalamus and peripheral clocks in peripheral tissues. According to the results of circadian transcriptomic studies in several tissues, the majority of rhythmic genes are expressed in a tissue-specific manner and are influenced by tissue-specific circadian rhythms. Here we review the diurnal variations of musculoskeletal functions and discuss the impact of the circadian clock on homeostasis in skeletal muscle and bone. Peripheral clocks are controlled by not only photic stimulation from the central clock in the SCN but also by external cues, such as feeding and exercise. In this review, we discuss the effects of feeding and exercise on the circadian clock and diurnal variation of musculoskeletal functions. We also discuss the therapeutic potential of chrono-nutrition and chrono-exercise on circadian disturbances and the failure of homeostasis in skeletal muscle and bone.

#### Edited by:

Etienne Challet, CNRS and University of Strasbourg, France

#### Reviewed by:

Stefano Schiaffino, Venetian Institute of Molecular Medicine, Italy Ke Ma, Baylor College of Medicine, USA

> \*Correspondence: Shigenobu Shibata<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Neuroendocrine Science, a section of the journal Frontiers in Neuroscience

Received: 28 November 2016 Accepted: 27 January 2017 Published: 14 February 2017

#### Citation:

Aoyama S and Shibata S (2017) The Role of Circadian Rhythms in Muscular and Osseous Physiology and Their Regulation by Nutrition and Exercise. Front. Neurosci. 11:63. doi: 10.3389/fnins.2017.00063 Keywords: circadian rhythm, clock gene, skeletal muscle, bone, chrono-exercise, chrono-nutrition

### INTRODUCTION

Various physiological functions, including the sleep wake cycle, body temperature, hormone secretion, and locomotor activity, exhibit circadian rhythms. This time-dependent regulation is driven by an internal circadian clock. In mammals, the circadian clock is divided into two parts, the master clock in the suprachiasmatic nucleus (SCN) of the hypothalamus and peripheral clocks in the peripheral tissues, such as the liver, skeletal muscle and so on, as well as brain areas other than the SCN. The master clock in the SCN acts as a time keeper in the whole body; thus, it integrates and entrains the peripheral circadian clocks by regulating neural and endocrine pathways, such as the sympathetic nervous system and glucocorticoid signaling (Schibler et al., 2003; Shibata, 2004). Light is the major entraining factor for the SCN. On the other hand, the peripheral clocks are entrained by not only the light-dependent regulation of the SCN but also scheduled feeding and scheduled exercise in an SCN-independent manner (Tahara and Shibata, 2013). The molecular mechanisms of circadian clock systems in mammals have been investigated since Clock (Circadian locomotor output cycles kaput) was discovered in 1997 (King et al., 1997). Circadian rhythm is produced by the transcriptional and translational regulation feedback loop of core clock genes, which include Bmal1 (Bain and muscle ARNT-like 1), Clock, Per1 (Period1), Per2, Cry1 (cryptochrome1), and Cry2. CLOCK and BMAL1 are transcriptional factors that have a basic helix-loop-helix PAS domain, and the heterodimer of CLOCK and BMAL1 proteins binds to an E-box binding element in the promoter regions of Per and Cry, and its binding activates the transcription of these genes (Gekakis et al., 1998). After transcription and translation, PER1/2 proteins localize to the cytoplasm and are phosphorylated by CKIε(Casein kinase Iε) (Lowrey et al., 2000). The phosphorylated PER1/2 proteins are not stable and are degraded by ubiquitination systems. CRY1/2 proteins in the cytoplasm are also degraded by ubiquitination systems via FBXL3 (F-box and leucine rich repeat protein 3) (Busino et al., 2007). The localization of CRY1/2 and PER1/2 to the cytoplasm promotes the formation of the PERs/CRYs/CKIε complex. This complex accumulates in the nucleus and suppresses the transcription of Pers and Crys by CLOCK and BMAL1. In addition, Rev-erbs (nuclear receptor subfamily 1, group D) and Rors (RAR-related orphan receptor) genes are the target genes of the complex of BMAL1 and CLOCK (Preitner et al., 2002; Sato et al., 2004). REV-ERBs and RORs act as the repressor and activator of Bmal1 and Clock transcription, respectively, via binding to an ROR-responsive element (Preitner et al., 2002; Sato et al., 2004). Recent circadian transcriptomic studies revealed that clock genes are expressed in several tissues, and that the expression of rhythmic genes in each tissue occurs in a tissue-specific manner (Miller et al., 2007; Zhang et al., 2014), suggesting that peripheral clocks in each tissue generate the biological rhythm of tissue-specific functions.

Clock genes are also expressed in skeletal muscle and bones of the skeleton, and it is thought that the clock genes regulate muscular- or osseous-specific biological functions (Dudek and Meng, 2014; Mayeuf-Louchart et al., 2015; Chatterjee and Ma, 2016; Yang and Meng, 2016). Skeletal muscle and bone have roles in not only the regulation of locomotion and postural support but also the control of nutritional homeostasis, such as maintaining glucose and calcium levels. Feeding and exercise stimulate these tissues and change their functions, including the maintenance of tissue mass and metabolism. Here, we review and discuss mainly two topics as follows: (1) the role of circadian rhythms in the biological functions of muscles and bone, and (2) the entrainment or regulation of the circadian clock or biological rhythm in skeletal muscle and bone by feeding and exercise.

### CIRCADIAN RHYTHM IN SKELETAL MUSCLES

Most mammalian cells express molecular clock genes and have a circadian clock system. Skeletal muscle cells also express molecular clock genes and show a circadian rhythm expression pattern. Using DNA microarray and RNA-seq methods Zhang et al. reported on the circadian transcriptome of 12 mouse organs and found that most rhythmic genes show an organ-specific pattern (Zhang et al., 2014). The circadian transcriptome of adult mouse skeletal muscle has also been identified by Miller et al., McCarthy et al., and Dyar et al. (McCarthy et al., 2007; Miller et al., 2007; Dyar et al., 2015). A small number, 3.4%, of the expressed genes in skeletal muscle show a circadian rhythm expression pattern (Miller et al., 2007). In addition, it has been reported that the number of rhythmic genes in skeletal muscle depends on the muscle fiber type; 684 rhythmic genes were found in the fast-twitch tibialis anterior muscle, while 1359 were identified in the slow-twitch soleus muscle (Dyar et al., 2015). The phase of expression of many rhythmic genes in skeletal muscle occurs at the mid-point of the subjective active phase (McCarthy et al., 2007; Miller et al., 2007). In liver, the phase expression of many rhythmic genes is different from that in skeletal muscle, the large cluster of phase expression occurs at the mid-point of subjective day and subjective night, respectively. One of the reasons for this may be that scheduled exercise can entrain the circadian clocks in skeletal muscles (Wolff and Esser, 2012), i.e., the phase of rhythmic gene expression in skeletal muscle may be regulated by the rhythm of locomotor activity. These results from circadian transcriptomic studies in skeletal muscle suggest circadian regulation of several muscle functions. For example, Myod1 (myogenic differentiation 1), Ucp3 (uncoupling protein 3), Atrogin1 (F-box protein 32), Myh1 (myosin heavy chain 1) are muscle-specific genes that play roles in myogenesis, muscle lipid utilization, protein metabolism and the organization of myofilaments, respectively. They have shown the circadian rhythms of gene expression, some of which have been shown to be directly regulated by clock genes (Andrews et al., 2010; Zhang et al., 2012). On the other hand, since the rhythmic expressions of Myod1 and Atrogin1 in ad libitum feeding mice are not observed in fasting mice (Shavlakadze et al., 2013), feeding and fasting patterns may be important factors for the rhythmic expression of muscle specific genes.

#### THE ROLE OF THE CIRCADIAN CLOCK IN MUSCLE MASS, MUSCLE STRENGTH, AND MYOFIBER TYPE

The roles of molecular circadian clocks in skeletal muscle mass, strength, and myofiber type are demonstrated with the use of molecular clock gene-deficient or mutant mice (**Table 1**). Kondratov et al. were the first to report on the effect of the molecular clock on muscle phenotype (Kondratov et al., 2006).   Yang et al. reported that body weight loss and early aging were not observed in tamoxifen-inducible Bmal1 knockout mice under the treatment of tamoxifen after muscle development (>3 months) (Yang et al., 2016). In addition, Dyar et al. showed that reductions in muscle weight and strength were not observed in inducible muscle-specific Bmal1 knockout mice once Bmal1 knockout was induced after development (Dyar et al., 2014). These reports suggest that Bmal1 expression during development is important for the maintenance of weight gain and muscular strength. In addition to the timing of Bmal1 expression, the role of intrinsic muscular Bmal1 in

#### TABLE 1 | Summary of the muscular and osseous phenotypes in clock gene mutant mice.

(Continued)

the muscle weight and strength has also been reported (Dyar et al., 2014; Schroder et al., 2015). Muscle-specific Bmal1 knockout mice showed an increase in muscle weight and a decrease in muscle strength, although the whole body Bmal1 knockout mice exhibited muscle loss. These reports indicate that while a circadian clock during development regulates muscle mass, an intrinsic skeletal muscular clock may not directly regulate it. Lack of muscle use and the lowering of daily locomotor activity levels reduce muscle mass (Powers et al., 2005). Disruption of Bmal1 in the whole body, and not in only skeletal muscle, induces the disappearance of an activity rhythm and a reduction in daily locomotor activity levels (Kondratov et al., 2006). In another report, musclespecific over expression of Bmal1 partially restores the Bmal1 knockout-induced reduction of activity levels and body weight loss without improving arrhythmic behavior (McDearmon et al., 2006). In addition, the locomotor activity of muscle-specific Bmal1 knockout mice, which show an increase in muscle mass, have a normal circadian rhythm, and is increased activity during the active phase (Dyar et al., 2014). These reports suggest that circadian clock-regulated activity levels may have an important role in the growth of skeletal muscle. Taken together, the results suggests that early aging and muscle weight loss in whole body Bmal1 knockout mice are controlled in non-myofiber cells by Bmal1 during development. Although the role of Bmal1 for the developing fetus has not been described, it is possible that the rhythmic expression-independent functions of Bmal1 in skeletal muscle growth are observed because Bmal1, just as other core clock genes, does not show rhythmic expression in the embryo (Dolatshad et al., 2010). Indeed, Lipton et al. found that the phosphorylation of BMAL1 by S6K1 (ribosomal S6 protein kinase 1) regulates translation, independently of its regulatory role in transcription (Lipton et al., 2015), suggesting that the phosphorylation of BMAL1 may affect the synthesis of muscle proteins, which is an important process for muscle growth. Further investigation is required of the role of molecular clocks in skeletal muscle growth.

Bmal1 knockout mice and Clock mutant mice have disrupted myofilament architecture in their skeletal muscle and decreased muscle strength at a single-fiber level (Andrews et al., 2010). The disruption of the myofilament architecture is observed in Myod knockout mice as well as Bmal1 knockout mice (Andrews et al., 2010). There is a binding site for the heterodimer of CLOCK and BMAL1 in the promoter region of Myod1 and this heterodimer directly regulates the rhythmic expression of Myod1, which is not seen in Bmal1 knockout mice or Clock mutant mice (Andrews et al., 2010; Zhang et al., 2012). MyoD plays a role in myogenesis, which includes the formation of myofibers from satellite cells and myoblast cells. Cahtterjee et al. demonstrated that BMAL1 is one of the key players of myogenesis (Chatterjee et al., 2013). Deficiency of Bmal1 in myoblast cells suppresses myogenesis related-gene expression, including Myod, Myf5 (myogenic factor 5), and Myogenin expression, and impairs the differentiation of myoblasts to myofibers (Chatterjee et al., 2013). Moreover, they have shown a role for the Wnt (wingless-type MMTV integration site family) pathway in Bmal1-induced myoblast differentiation, since the circadian regulation of Wnt pathwayrelated genes is regulated by Bmal1 (Chatterjee et al., 2013). In addition, they have demonstrated that Bmal1 in skeletal muscle also promotes skeletal muscle regeneration via satellite cell proliferation using in vivo muscle injury models (Chatterjee et al., 2015). However, Schiaffino et al. reported that Myod1 expression in the skeletal muscle of muscle-specific Bmal1 knockout mice showed rhythmicity and was increased during the active phase (Schiaffino et al., 2016), while a previous study showed that whole body Clock∆<sup>19</sup> mice did not show rhythmicity (Andrews et al., 2010). As one of the reasons, the up regulation of Myod1 in muscle-specific Bmal1 knockout mice may be due to the increase in locomotor activity during the active phase (Dyar et al., 2014); this hypothesis is supported by the finding that mechanical stimulation, such as during exercise, increases Myod expression (Legerlotz and Smith, 2008). Indeed, a stimulation such as fasting, which might cause hypolocomotion, down regulated the rhythmic expression of Myod1 under conditions when the rhythmic expression of Bmal1 is maintained (Shavlakadze et al., 2013). As another reason, the rhythmic expression of Myod1 in muscle-specific Bmal1 knockout mice may be because of the tissue-specific promoter driving genetic ablation. Myosin light chain 1f (Mlc1f) promoter was used to drive the ablation of muscle-specific Bmal1 (Dyar et al., 2014). Mlc1f promoter is active in mature myocytes, but not in muscle precursor cells, such as satellite cells (Lee et al., 2012), suggesting that deficiency of Bmal1 occurs only in mature myocytes. Myod is expressed not only in muscle precursor cells, such as activated satellite cells and myoblasts, but also in the mature muscle myocytes, albeit at low levels (Voytik et al., 1993; Hughes et al., 1997). The transcriptional regulation of Myod occurs through two elements, the "core enhancer" and a bipartite element containing the "distal regulatory region" and the "proximal regulatory region" (Charge et al., 2008). Myod expression in myoblasts is activated through the core enhancer, and its expression in mature myocytes is up regulated through the bipartite element (Charge et al., 2008). In addition, BMAL1 and CLOCK bind to the core enhancer, but not the bipartite element (Andrews et al., 2010). The rhythmic expression of Myod1 in muscle-specific Bmal1 knockout mice may reflect its rhythmic expression in activated satellite cells and myoblasts, which are included in adult skeletal muscles. Therefore, these reports indicate that clocks in myoblasts and satellite cells and locomotor activity rhythms regulate skeletal muscle functions, including the formation and growth of skeletal muscle, via skeletal muscle specific regulators such as MyoD.

Skeletal muscle fiber is divided into the two types, slow-twitch fiber and fast-twitch fiber, according to the myosin heavy chain isoform composition (Schiaffino and Reggiani, 1996). The slowtwitch fiber is mainly composed of myosin heavy chain isoforms I and IIa. It has high oxidative capacity and mitochondrial volume. On the other hand, the fast-twitch fiber is mainly composed of myosin heavy chain isoforms IIx and IIb. Some circadian clock genes regulate fiber type composition in skeletal muscle. Muscle-specific Bmal1 knockout mice show a shift in their fiber type from slow- to fast-twitch compared with wild type mice (Dyar et al., 2014; Schroder et al., 2015). Deficiency of Reverbα, a repressor of Bmal1 transcription, induces a shift in fiber type from fast- to slow-twitch in the slow-twitch soleus muscle, however, its effects are small (Pircher et al., 2005). On the other hand, Woldt et al. showed that fiber type did not shift significantly in the wild type and Rev-erbα mice, although the genetic expression of slow-twitch fiber markers was lower in muscle from Rev-erbα mice than in wild type mice (Woldt et al., 2013). The reason for this discrepancy may be in the independent line of the global genetic targeting of the Rev-erbα allele. While there are few reports regarding the relationship between molecular clocks and muscle fiber types, Bmal1 has the potential to shift fiber types from fast- to slow-twitch. In contrast, Rev-erbα, the repressor of Bmal1, has opposing effects on fiber types.

### THE ROLE OF THE CIRCADIAN CLOCK ON LIPID AND CARBOHYDRATE METABOLISM IN SKELETAL MUSCLE

As mentioned previously, there is a lot of evidence regarding the regulatory effects of circadian rhythm on lipid and carbohydrate metabolism in the liver (Tahara and Shibata, 2016). In skeletal muscle, circadian rhythm-regulated lipid and carbohydrate metabolism may be due to intrinsic molecular clocks or could be a response to a behavior such as feeding/fasting or neural and hormonal cues. Hodge et al. have reported that Bmal1 in skeletal muscle-regulated genes is involved in the utilization and storage of energy substrates, independent of circadian activity (Hodge et al., 2015). For example, the peak time of lipogenic and lipolytic gene expression occurs at the end of the active phase and at the middle of the inactive phase, respectively. Carbohydrate catabolism and storage peak at the beginning of the active phase and in the middle of the active phase, respectively. These substrate-dependent peak times for metabolic genes in skeletal muscle are partially controlled by an intrinsic skeletal muscle molecular clock (Hodge et al., 2015). Other reports have demonstrated that insulin-stimulated glucose uptake, its metabolism, and its related factors are down-regulated in skeletal muscle-specific Bmal1 knockout mice (Dyar et al., 2014; Harfmann et al., 2016). Glucose transporter 4 (GLUT4) is a key molecule for glucose transport in skeletal muscle. Insulin translocates GLUT4 from the cytoplasm to the plasma membrane via the activation of the insulin signaling pathway, which includes TBC1 domain family, member 1 (TBC1D1). Quantities of GLUT4 and TBC1D1 show a diurnal change, with increased levels during the active phase and decreased levels during the inactive phase (Dyar et al., 2014). These diurnal changes are not observed in muscle-specific Bmal1 knockout mice (Dyar et al., 2014) and a deficiency in Bmal1 reduces levels of these molecules throughout the day (Dyar et al., 2014; Harfmann et al., 2016). On the other hand, the deficiency of muscular Bmal1 did not affect insulin signaling, including the phosphorylation of Akt (Dyar et al., 2014). In addition to affecting the transport of glucose, a deficiency in Bmal1 causes dysregulation of glycolysis and glucose oxidation via the inactivity of metabolic enzymes, such as Hexokinase 2 (HK2) and Pyruvate dehydrogenase (PDH), suggesting abnormal glucose metabolism (Dyar et al., 2014; Harfmann et al., 2016). Owing to the dysregulation of the glycolytic pathway, a deficiency of Bmal1 in skeletal muscle increases levels of the metabolites related to the pentose phosphate pathway, the polyol pathway, and glucuronic acid pathway (Dyar et al., 2014). These results, which relate not only to the transcriptome, but also to protein levels and metabolite levels, strongly support the notion that a major physiological role of the muscle clock is to prepare for the transition from the rest/fasting phase to the active/feeding phase, when glucose becomes the predominant fuel for skeletal muscle. Another investigation, using C2C12 myotubes, showed the role of muscle clocks in insulin sensitivity through Sirtuin 1 (Sirt1) (Liu et al., 2016). In this report, the knockdown of Clock and Bmal1 caused insulin resistance via Sirt1 (Liu et al., 2016). The

muscular clock-regulated Rev-erbα also shows high expression levels in slow-twitch fiber type muscle, such as soleus muscle, and regulates lipid uptake and oxidative capacity in skeletal muscle by controlling mitochondrial biogenesis and autophagy (Woldt et al., 2013). A deficiency in Rev-erbα reduces oxidative capacity, which causes exercise intolerance (Woldt et al., 2013). RORα is highly expressed in skeletal muscles (Becker-Andre et al., 1993). A deficiency in Rorα inhibits the expression of genes involved in lipid homeostasis in skeletal muscle cells (Lau et al., 2004). In particular, RORα directly regulates the transcription of Carnitine palmitoyltransferase-1 (Cpt-1) and Caveolin-3 (Cav3) (Lau et al., 2004). These data suggest that the molecular clock generates the circadian rhythm of metabolism in skeletal muscle and that the disruption of circadian rhythm occurs owing to metabolic dysfunction in skeletal muscle.

### THE PHASE-SHIFTED EFFECTS OF EXERCISE AND FEEDING TIME ON MUSCULAR CIRCADIAN CLOCKS

The hierarchy of tissue-specific clocks exists in the mammalian circadian clock system.  Brain and peripheral clocks are under the control of the master clock. Classically, the master clock is entrained by a photic cue, such as the light-dark cycle, and it is thought that peripheral clocks are also regulated by the central clock in the SCN. On the other hand, peripheral clocks are controlled in a central-clock-independent manner under certain conditions, such as scheduled exercise and restricted feeding (Damiola et al., 2000; Stokkan et al., 2001). Exercise is one of the non-photic phase-shifting cues (**Figure 1**). Several studies have shown scheduled exercise during the daytime or subjective daytime advancing of the phase of circadian rhythms in rodents (Reebs and Mrosovsky, 1989; Marchant and Mistlberger, 1996). Wolff et al. have shown the phase advanced effects of scheduled exercise in ex-vivo experiments (Wolff and Esser, 2012). After scheduled voluntary or involuntary (treadmill running) exercise during the daytime for 4 weeks, the phase expression of Per2::Luc advances in skeletal muscle and lung but not in the SCN (Wolff and Esser, 2012). Scheduled exercise also accelerates reentrainment in mouse skeletal muscle and lung but not liver or the SCN to a new light-dark cycle (Yamanaka et al., 2008). In addition, Yamanaka et al. have shown that exercise-induced re-entrainment depends on the timing of exercise and on the peripheral tissues (Yamanaka et al., 2016). In the phase advanced light/dark cycle condition, wheel running at the beginning of the active phase (onset) accelerates the re-entrainment of the skeletal muscle clock but not the SCN. On the other hand, in the phase delayed light/dark cycle condition, wheel running at the end of the active-phase (offset) interferes with the re-entrainment of the skeletal muscle clock. These reports have indicated the potential of exercise to induce a phase shifting effect in skeletal muscle with the use of ex vivo experiments. In recent years, we have reported that exercise advances the phase of circadian rhythm in peripheral clocks, such as liver and gastrocnemius muscle, in in vivo experiments (Sasaki et al., 2016). In this report, the exercise-induced phase advance of the liver clock was stronger following forced exercise rather than voluntary exercise, although the difference between the two kinds of exercise to cause entrainment was not observed in the skeletal muscle clock. While the mechanism of exercise-induced phase shifting has not been fully elucidated, the role of some factors has been proposed to date.  In fact, we have shown that forced treadmill running increases serum corticosterone and tissue norepinephrine levels and their elevations play an important role in the forced exercise-induced phase shifting of peripheral clocks (Sasaki et al., 2016). The effects of inactivity on peripheral clocks have also been reported (Dyar et al., 2015; Nakao et al., 2015). Nakao et al. showed that the expression of Bmal1, Per1, Per2, Rora, Nr1d1, and Dbp were decreased, and the expression of Clock was increased, in inactive muscles by denervation of the sciatic nerve (Nakao et al., 2015). In denervated muscles, Bmal1 and Dbp showed phase advance, compared with the contralateral muscle. The denervation-induced phase advance of muscular clocks has also reported by Dyar et al. (2015). Although the loss of muscle activity did not completely abolish the rhythms of muscular clocks, the rhythmic expressions of many cyclic genes were altered by denervation (Dyar et al., 2015; Nakao et al., 2015). These data suggest that physical activity affects muscular clocks. Further investigation is needed into the role of neuronal signals and physical activity in the direct and indirect regulation of muscular clocks; physical activity causes several physiological changes, such as to body temperature and hormonal status, which are known to affect the peripheral clocks (Tahara et al., 2017). The phase shifting effects of exercise have been observed in human studies (Van Reeth et al., 1994; Buxton et al., 1997, 2003; Miyazaki et al., 2001; Barger et al., 2004). For example, Yamanaka et al. reported that the sleep-wake cycle but not melatonin rhythm is accelerated by exercise under dim light conditions and a phaseadvanced sleep schedule (Yamanaka et al., 2010). Zambon et al. reported that resistant exercise changed the gene expression of the circadian clock in the skeletal muscle of humans (Zambon et al., 2003). These reports indicate that exercise acts as a potential entrainer of skeletal muscle clocks in humans as well as rodents.

As mentioned before, the peripheral clock is regulated by restricted feeding (Damiola et al., 2000; Stokkan et al., 2001). To date, restricted feeding shifts and entrains the liver clock in an SCN-independent manner (Hara et al., 2001; Tahara and Shibata, 2014). It is thought that insulin is the hormone responsible for the feeding-induced phase advance of the peripheral clock (Kuriyama et al., 2004; Tahara et al., 2011; Dang et al., 2016). Skeletal muscle is one of the major insulin sensitive organs and it controls whole body blood glucose levels via insulinstimulated glucose uptake (DeFronzo, 1988). This evidence

gives the expectation that clocks in insulin-sensitive organs are regulated by restricted feeding, in other words restricted feeding regulates the peripheral clock not only in the liver but also in skeletal muscle. However, Guo et al. reported tissue-dependent regulation of the peripheral clock by hormonal cues in SCNlesioned mice (Guo et al., 2005). When circulating factors in SCN-lesioned mice and intact mice are shared by parabiosis, the arrhythmic clock gene expression in the liver of SCN-lesioned mice is restored. On the other hand, the disturbance of the skeletal muscle clock is not restored by parabiosis, indicating that circulating factors, such as hormones, are not important for the regulation of the skeletal muscle clock. In fact, some reports have shown that the expression of core clock genes, such as Bmal1 and Per2, is synchronized by day time feeding in the liver but not skeletal muscle (Reznick et al., 2013; Yasumoto et al., 2016). On the other hand, some rhythmic genes, such as Pdk4 (pyruvate dehydrogenase kinase 4) and Ucp3, in skeletal muscles show arrhythmic or dampened expression patterns owing to day time feeding. These data suggest that restricted feeding does not have a major effect on the skeletal muscle clock. However, further evidence is needed because reports about the effect of restricted feeding on circadian rhythm in skeletal muscle have been limited.

### CIRCADIAN RHYTHM OF BONE FUNCTIONS

Bone is one of the organs to play a major role in the storage of calcium and phosphorus. Three kinds of cells exist in bone, namely osteocytes, osteoblasts, and osteoclasts, (Prideaux et al., 2016). More than 90% of cells in bone are osteocytes, and these cells play a role in the storage of bone matrix (Prideaux et al., 2016), whereas osteoblasts and osteoclasts regulate bone remodeling. Bone homeostasis is regulated by the balance between the production of bone matrix by osteoblasts, called bone formation, and the breakdown of bone matrix by osteoclasts, called bone resorption (Rodan and Martin, 2000). Not only the process of bone remodeling but also serum concentrations of some hormones regulating bone metabolism show diurnal variation (Dudek and Meng, 2014). For example, it has been reported that serum concentrations of calcitonin, calcium, osteocalcin, parathyroid hormone Ctelopeptide, skeletal alkaline phosphatase, and tartate-resistant acid phosphatase show diurnal variation (Greenspan et al., 1997; Srivastava et al., 2001; Shao et al., 2003; Yang and Meng, 2016). In addition, Zvonic et al. reported on the circadian transcriptome in calvarial bone (Zvonic et al., 2007). In this report, more than 26% of genes expressed in calvarial bone show a rhythmic expression (Zvonic et al., 2007). Interestingly, among the rhythmic genes in calvarial bone 64% of genes do not exhibit rhythmic expression in liver, brown adipose tissue or white adipose tissue. In other words, a lot of rhythmic gene expression is controlled in a tissue-specific manner. For example, calcium channels, ADAMs (a disintegrin and metalloproteinases), FGFs (fibroblast growth factors), and Runxs (runt related transcription factors), which have critical roles in bone formation, bone remodeling, bone metabolism regulated cytokines and bone-specific transcriptional factors, show diurnal gene expression patterns (Zvonic et al., 2007). The transcriptomic results from calvarial bone suggest that the process of metabolism in several bones is regulated by circadian rhythm.

Indeed, it has been reported that some molecular clocks regulate and maintain bone homeostasis (**Table 1**). Fu et al. demonstrated that both Per1 deficient and Per2 PAS domain mutant mice have high bone mass (Fu et al., 2005). This phenotype is also exhibited by other clock gene modified mice, i.e., Per1 and 2 double knockout mice, Cry1 and 2 double knockout mice, Bmal1 knockout mice, Per2 mutant mice and Cry2 knockout mice (Fu et al., 2005; Maronde et al., 2010). These reports have suggested that the molecular clock is a negative regulator of bone mass. The high bone mass induced by molecular clock-deficiency is controlled by osteoblast differentiation via leptin-Adrβ2 (adrenoceptor beta 2)-Creb1 cascades (Fu et al., 2005). Indeed, leptin levels are elevated in mice with both genotypes, Per1−/<sup>−</sup> and Per2−/<sup>−</sup> mutant mice, and the elevation of leptin levels may be involved in osteoblast differentiation via Adrβ2 activation. In fact, it has been reported that both leptin deficient mice and Adrb2 deficient mice display a low bone mass phenotype (Ducy et al., 2000; Takeda et al., 2002). Therefore, clock deficiency-induced high bone mass is involved in leptin-dependent sympathetic activation via the Adrβ2. In contrast, Samsa et al. reported that deficiency of Bmal1 in mice results in age-related bone loss (Samsa et al., 2016), which is in contrast to the results of a previous report (Fu et al., 2005). The different bone phenotypes in these reports may be due to the age differences of the mice. Pers, Crys, or Bmal1 deficiencyinduced high bone mass is observed during young and adolescent ages (Fu et al., 2005; Maronde et al., 2010). On the other hand, Bmal1 deficiency-induced low bone mass is observed at geriatric age (Samsa et al., 2016). Bmal1 knockout mice also show ageassociated phenotypes such as body weight loss and sarcopenia (Kondratov et al., 2006). In addition, abnormal bone calcification and arthropathy in Bmal1 knockout mice did not replicate in tamoxifen-inducible Bmal1 knockout mice after development (Yang et al., 2016). While the reason for the opposing effects of Bmal1 on bone mass have not been well elucidated, the age-dependent effects of Bmal1 on bone mass may be one explanation.

The effects of circadian clock genes on bone resorption have been examined using osteoclast-specific Bmal1 knockout mice (Xu et al., 2016). The osteoclast-specific Bmal1 knockout mice have a high bone mass. The upregulation of Nfatc1 via the direct regulation of the heterodimer of BMAL1 and CLOCK induces osteoclast differentiation and then reduces bone mass (Xu et al., 2016). This report suggests that the osteoclast clock genes reduce bone mass via the activation of osteoclast differentiation.

Interestingly, in recent years it has been reported that skeletal muscle clocks are linked to the maintenance of bone homeostasis. Schroder et al. reported that muscle specific-Bmal1 knockout mice show thickening of the distal tibia (Schroder et al., 2015). Muscle and bone may communicate with each other via endocrine factors and mechanical loading (Karsenty and Olson, 2016). For example, myostatin, one of the hormones secreted from skeletal muscle, regulates osteoclast differentiation (Dankbar et al., 2015). While the mechanism of interaction of these tissues has not been fully elucidated yet, the skeletal muscle clock is closely linked to the maintenance of bone homeostasis.

### REGULATION OF THE OSSEOUS CIRCADIAN RHYTHM BY INTERNAL AND EXTERNAL CUES

The circadian rhythms of bone functions are controlled by internal or external cues to maintain a balance between bone formation and bone resorption. The master clock in the SCN controls the peripheral circadian clocks via output signals such as hormonal and sympathetic nervous system signaling. Leptin is known to be one of the entraining factors of the circadian clock in osteoblasts (Fu et al., 2005) (**Figure 1**). Leptin activates the sympathetic nervous system via the Adrβ2, of which Creb1 is a downstream factor (Fu et al., 2005). These reports indicate that the osteoblast clock is regulated by sympathetic nervous system signals. In addition, in human osteoblasts, treatment with a β-adrenergic receptor agonist, isoprenaline, or synthetic glucocorticoids, dexamethasone, induces the circadian expression of clock genes (Komoto et al., 2012), suggesting that both sympathetic nervous system and glucocorticoid signals are involved in the circadian rhythm in osteoblasts and act as mediators from the SCN to osteoblasts. On the other hand, in osteoclasts, glucocorticoid signals rather than sympathetic signals have the ability to regulate rhythmic gene expression. Fujihara et al. demonstrated that the rhythmic expression of clock genes is changed by stimulation with dexamethasone but not isoprenaline (Fujihara et al., 2014). Osteoclast-specific genes, Nfatc1 (nuclear factor of activated T cells 1) and Ctsk (cathepsin K), show rhythmic expressions (Fujihara et al., 2014). The rhythmic expression of Nfatc1 participates in Bmal1-regulated bone resorption, as described above (Xu et al., 2016). These rhythmic expressions are dampened in cancellous bone from adrenalectomized mice, and glucocorticoid injection restores the rhythmic expression of these genes (Fujihara et al., 2014). This evidence indicates that glucocorticoid signaling acts as a mediator between the SCN and osteoclasts for the synchronization of circadian rhythms in osteoclasts.

Feeding and fasting regulate the diurnal variation of the bone resorption marker serum C-telopeptide fragments of collagen type 1 degradation (s-CTx) (Gertz et al., 1998) (**Figure 1**). The levels of s-CTx in humans are higher during early morning, from 05:00 to 08:00, and lower in the late afternoon, from 12:00 to 16:00. The levels of s-CTx show diurnal variation under normal feeding conditions, such as the consumption of breakfast, lunch and dinner, while fasting dampens the diurnal variation of s-CTx (Bjarnason et al., 2002; Qvist et al., 2002). In addition, the feeding-induced generation of this diurnal variation is also observed in the intake of glucose, protein, and fat (Bjarnason et al., 2002). The diurnal variation of levels of s-CTx under normal feeding conditions only occurs during breakfast but not lunch or dinner (Bjarnason et al., 2002). A reason for this may be that fasting plays a role in the feeding-induced generation of the s-CTx rhythm, since the fasting period from breakfast to lunch and from lunch to dinner is shorter than that from dinner to breakfast. These reports suggest that the feeding and fasting rhythm is an important factor in regulating and maintaining the circadian rhythm of bone resorption, although the mechanism for the preventative effects of food intake on bone resorption has not been well investigated.

Mechanical loading, such as exercise, protects against agerelated bone loss, whereas unloading, such as bed rest, induces bone loss (Qi et al., 2016). Some reports have shown the effects of unloading on circadian rhythms of bone resorption markers (Halloran et al., 1988; Pedersen et al., 1995; Kim et al., 2000). Pedersen et al. showed that 5 days of bed rest in healthy women did not change the circadian rhythm of s-CTx or other serum bone resorption markers, such as alkaline phosphatase and osteocalcin (Pedersen et al., 1995). No effects of unloading on diurnal variation of bone resorption markers have been shown by other reports (Halloran et al., 1988; Kim et al., 2000). These reports suggest that a common physical activity rhythm, which also includes low intensity exercise, does not have the potential to regulate the circadian rhythm of bone resorption. However, moderate or high intensity scheduled exercise may affect the circadian variation of bone metabolism, since mechanical loading influences bone resorption and bone formation, as well as circulating factors such as glucocorticoids and parathyroid hormone (Fragala et al., 2011; Gardinier et al., 2015; Qi et al., 2016; Sasaki et al., 2016).

### PERSPECTIVES

Based on recent findings, circadian rhythms in skeletal muscle and bone maintain their homeostasis. The disruption of muscle clocks occurs owing to dysregulation of whole body glucose metabolism (Harfmann et al., 2016) and bone clocks could be negative regulators of bone mass through the inhibition of bone formation and the activation of bone resorption (Fu et al., 2005; Xu et al., 2016). These finding suggest that disturbances of circadian rhythms by social or environmental factors, such as shift work, may result in dysfunctions of skeletal muscle and bone. In epidemiological studies, the prevalence of metabolic syndrome, osteoporosis and bone fractures is increased in shift workers (Feskanich et al., 2009; Pietroiusti et al., 2010; Quevedo and Zuniga, 2010). In addition, it has been reported that long term constant light exposure reduces muscle strength and bone mass in mice (Lucassen et al., 2016). These findings indicate that the regulation of circadian rhythms in skeletal muscle and bone by external cues, such as feeding and exercise, are important for the maintenance of homeostasis in these tissues, since circadian rhythm in these tissues can be entrained or regulated by the feeding/fasting rhythm and the physical activity rhythm, including scheduled exercise. We report that disturbances of circadian rhythms in peripheral clocks under constant light conditions are partially improved by scheduled feeding and scheduled exercise, although the circadian clocks in skeletal muscle and bone have not been measured (Hamaguchi et al., 2015). In addition, the reductions in both mitochondrial content and exercise tolerance in the skeletal muscle of Clock mutant mice are restored by endurance exercise

#### REFERENCES

training (Pastore and Hood, 2013), suggesting that exercise can prevent circadian disturbance-induced muscular dysfunctions. In addition, some reports have shown the importance of the circadian timing of exercise for the prevention of diseases. Schroeder et al. have reported that rhythmic deficits observed in vasointestinal polypeptide-deficient mice are improved by wheel running exercise (Schroeder et al., 2012). Interestingly, greater preventative effects are observed when wheel running occurs at the end of the active phase but not at the beginning of the active phase (Schroeder et al., 2012). We have also reported that wheel running exercise at the end of the active phase has more preventative effects on high fat diet-induced obesity than that occurring at the beginning of the active phase (Sasaki et al., 2014). In addition, exercise in the morning but not in the afternoon or evening increases fat oxidation over 24 h in healthy humans (Iwayama et al., 2015). These reports suggest the importance of scheduled feeding or exercise for skeletal muscle and bone health. However, the importance of the circadian timing of exercise and nutritional intake for muscular and osseous health has not been well elucidated. Further advanced evidence is required and it is expected to lead to a better understanding of the mutual interaction between the circadian clock and muscle/bone.

#### AUTHOR CONTRIBUTIONS

SA and SS were involved in conceptualizing and writing the manuscript.

#### FUNDING

This work was supported by the council for Science, Technology, and Innovation, SIP, "Technologies for creating next-generation agriculture, forestry, and fisheries" (funding agency: Bio-oriented Technology Research Advancement Institution, NARO) (SS).

in response to S6K1-mediated phosphorylation. Cell 161, 1138–1151. doi: 10.1016/j.cell.2015.04.002

sleep-wake cycle but not of plasma melatonin rhythm to 8-h phase-advanced sleep schedule. Am. J. Physiol. Regul. Integr. Comp. Physiol. 298, R681–R691. doi: 10.1152/ajpregu.00345.2009

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Aoyama and Shibata. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# High-Resolution Recording of the Circadian Oscillator in Primary Mouse **α**- and **β**-Cell Culture

*Volodymyr Petrenko1,2,3,4, Yvan Gosmain4,5 and Charna Dibner 1,2,3,4\**

*1Endocrinology, Diabetes, Hypertension and Nutrition Division, Department of Specialties of Medicine, University Hospital of Geneva, Geneva, Switzerland, 2Faculty of Medicine, Department of Cell Physiology and Metabolism, University of Geneva, Geneva, Switzerland, 3Institute of Genetics and Genomics in Geneva (iGE3), Geneva, Switzerland, 4Diabetes Center of the Faculty of Medicine, University of Geneva, Geneva, Switzerland, 5Molecular Diabetes Laboratory, Endocrinology, Diabetes, Hypertension and Nutrition Division, Faculty of Medicine, Department of Specialties of Medicine, University Hospital of Geneva, University of Geneva, Geneva, Switzerland*

Circadian clocks have been developed in evolution as an anticipatory mechanism allowing for adaptation to the constantly changing light environment due to rotation of the Earth. This mechanism is functional in all light-sensitive organisms. There is a considerable body of evidence on the tight connection between the circadian clock and most aspects of physiology and metabolism. Clocks, operative in the pancreatic islets, have caught particular attention in the last years due to recent reports on their critical roles in regulation of insulin secretion and etiology of type 2 diabetes. While β-cell clocks have been extensively studied during the last years, α-cell clocks and their role in islet function and orchestration of glucose metabolism stayed unexplored, largely due to the difficulty to isolate α-cells, which represents a considerable technical challenge. Here, we provide a detailed description of an experimental approach for the isolation of separate mouse α- and β-cell population, culture of isolated primary α- and β-cells, and their subsequent long-term high-resolution circadian bioluminescence recording. For this purpose, a triple reporter *ProGlucagon*-*Venus*/RIP-*Cherry*/*Per2*:*Luciferase* mouse line was established, carrying specific fluorescent reporters for α- and β-cells, and luciferase reporter for monitoring the molecular clockwork. Flow cytometry fluorescence-activated cell sorting allowed separating pure α- and β-cell populations from isolated islets. Experimental conditions, developed by us for the culture of functional primary mouse α- and β-cells for at least 10 days, will be highlighted. Importantly, temporal analysis of freshly isolated α- and β-cells around-the-clock revealed preserved rhythmicity of core clock genes expression. Finally, we describe the setting to assess circadian rhythm in cultured α- and β-cells synchronized *in vitro*. The here-described methodology allows to analyze the functional properties of primary α- and β-cells under physiological or pathophysiological conditions and to assess the islet cellular clock properties.

#### *Edited by:*

*Andries Kalsbeek, Academic Medical Center, Netherlands*

#### *Reviewed by:*

*Aleksey V. Matveyenko, Mayo Clinic Minnesota, USA Elaine Vieira, Researcher at Universidade católica de Brasília, Brazil*

> *\*Correspondence: Charna Dibner<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Neuroendocrine Science, a section of the journal Frontiers in Endocrinology*

*Received: 20 January 2017 Accepted: 24 March 2017 Published: 07 April 2017*

#### *Citation:*

*Petrenko V, Gosmain Y and Dibner C (2017) High-Resolution Recording of the Circadian Oscillator in Primary Mouse α- and β-Cell Culture. Front. Endocrinol. 8:68. doi: 10.3389/fendo.2017.00068*

Keywords: mouse pancreatic islet, **α**- and **β**-cells, primary culture, *in vitro* synchronization, circadian bioluminescence

## INTRODUCTION

The circadian system represents a complex anticipatory mechanism developed during evolution in nearly all organisms, allowing to coordinate a plethora of physiological functions to the daily changes of geophysical time. Within this system, a master pacemaker in the hypothalamus orchestrates subsidiary oscillators situated in peripheral organs (1). In fact, myriads of these self-sustained and cell-autonomous oscillators are operative in most cells of the body (2, 3). The molecular composition of central and peripheral oscillators is identical, and it relies on primary and secondary feedback loops of transcription and translation of key core clock components (4). The primary loop comprises the positive limb transcription factors CLOCK and BMAL1, which induce expression of the negative limb elements PERIODS and CRYPTOCHROMES (5). Recent studies provide increasing evidence for a tight connection between the circadian system and metabolism, linking metabolic diseases to circadian misalignments associated with modern life-style, including frequent jetlag, shifted work schedules, and chronic social jetlag (4, 6–10). Studies in clock-deficient genetic rodent models suggest that a number of metabolic defects develop in mice that are deficient for one or two core clock components (11, 12). For instance, *Clock* mutant mice develop hyperphagia, obesity, hyperglycemia, and hypoinsulinemia (12).

There is an increasing evidence for the essential roles of the peripheral circadian clocks operative in endocrine tissues for their transcriptional and functional regulation (13–15). Indeed, most of the hormones, including myokines and adipokines, are secreted in a circadian manner and regulated by respective cell-autonomous oscillators (16, 17). Such cell-autonomous clocks have been recently characterized in pancreatic islets in mice (11, 18) and in humans (18–20). Loss of islet clock function in islet-specific *Bmal1* KO mouse models, either induced during development or in the adult age, resulted in the early onset of type 2 diabetes (T2D) in these mice (11, 18, 21). Moreover, siClock-mediated clock perturbation in adult human islet cells caused disruption in basal and glucose induced insulin secretion by these cells *in vitro* (20). Taken together, these data suggest that circadian oscillators operative in islet cells play an important role in regulating these cell function.

So far, most of the research works were conducted on whole islets, or on insulin secreting β-cells, representing about 80% of total islet cells in mice (22). Therefore, the circadian physiology of glucagon secreting α-cells stayed largely unexplored, due to the difficulty to identify these cells within the complex three dimensional islet structure and to isolate them due to their low abundance (less than 20% of the mouse islet cell population). In an attempt to fill this gap, we hereby report an experimental approach, which allows to (1) efficiently isolate nearly pure populations of mouse α- and β-cells; (2) establish and maintain mouse α- and β-cell primary cultures; (3) study endocrine function of separated α- and β-cells; and (4) assess the circadian properties of primary α- and β-cells, utilizing high-resolution circadian bioluminescence monitoring in living cells synchronized *in vitro*.

#### MATERIALS AND METHODS

#### Animal Care and Reporter Mouse Strain

For all experiments a triple reporter mouse strain *ProGlucagon*-*Venus*/RIP-*Cherry*/*Per2*:*Luciferase* (*ProGcg*-*Venus*/RIP-*Cherry*/ *Per2*:*Luc*) was derived by crossing the *ProGlucagon*-*Venus (ProGcg-Venus)* reporter mouse (23) with Rat *Insulin2* promoter (RIP)-*Cherry* (RIP-*Cherry*) (24) and *Period2:Luciferase*  (*Per2:Luc*) mice (25). *ProGcg-Venus* and *RIP-Cherry* reporters exhibit a high specificity for α- and β-cells, respectively, while the fusion protein PER2:Luciferase, encoded by *Per2:Luc*, is a circadian reporter functionally indistinguishable from the wildtype PER2 protein. The overview of the experimental procedures is illustrated in **Figure 1**. All experiments were conducted on male mice aged 7–16 weeks under standard animal housing conditions comprising *ad libitum* access to food and water and 12 h light/12 h dark cycles. Islet isolations were performed during morning hours (07:00 a.m.–12:00 a.m.). To study circadian rhythms in freshly isolated α- and β-cells *in vivo*, mice were subjected to night-restricted feeding 2 weeks prior to the experiment and during the entire period of sample collection as described previously (26), with half of the animals entrained with inversed light-dark and feeding cycles during the same period.

#### Pancreatic Islet Isolation and Separation of **α**- and **β**-Cells

Islets of Langerhans were isolated by a standard procedure based on collagenase (Type XI, Sigma) digestion of pancreas followed

secretion]. Scale bar = 200 µm.

by Ficoll purification (27). Briefly, 2 ml of a 2-mg/ml collagenase solution [diluted in Hanks Balanced Salt Solution (HBSS) (Gibco), 4.2 nM sodium bicarbonate, 1 M CaCl2, and 1 M HEPES (Gibco)] was injected retrogradely through the ampulla of Vater into the exocrine part of the pancreas, which was dissected and digested in a water bath at 37°C (11–14 min, exact time needs to be calibrated for every new batch of enzyme). The digested mixture was then washed in HBSS, supplemented with 0.3% free fatty acid bovine serum albumin (BSA) (Sigma) with subsequent Ficoll gradient centrifugation (densities: 1.108, 1.096, and 1.069). Islets, collected from the gradient, were then washed in HBSS supplemented with 0.3% BSA. Islet cells were gently dissociated by trypsin (Gibco), re-suspended in KRB solution (pH 7.4, BSA, 1.4 mM glucose and 0.5 mM EDTA), and filtered through a 70-µm cell strainer (Falcon). α- and β-cell populations were separated by flow cytometry fluorescence-activated cell sorting [FACS; Astrios Sorter (Beckman Coulter)], based on fluorescence intensity and wavelength, single cell nature, size, and viability. Purity of the sorted cells was further assessed by a second, additional FACS analysis. The viability of treated cells was evaluated by FACS with DRAQ7™ (Abcam) or DAPI (ThermoFisher) dyes.

#### Immunohistochemistry

Sorted cells were plated on 35 mm dishes and fixed in 4% paraformaldehyde during 30 min. For immunohistochemistry, the fixed cells were incubated with mouse anti-Glucagon (1:500, Sigma) or guinea-pig anti-Insulin (1:500, ThermoFisher Scientific) primary antibodies. The signal was revealed by Alexa-Fluor anti-mouse 568 or Alexa-Fluor anti Guinea-pig 488 secondary antibodies, respectively (1:1,000, Molecular Probes). Bright-field and fluorescent images were obtained with EvosTL fluorescent microscope (ThermoFisher) using 4×, 10×, or 40× objectives.

#### *In Vitro* Islets/Islet Cell Culture

For the *in vitro* culture, intact islets or sorted cells were recovered in RPMI 1640 complete medium (11.2 mM glucose, 110 µg/ml sodium pyruvate) supplemented with 10% fetal calf serum, 110 U/mL penicillin, 110 µg/ml streptomycin, and 50 µg/ml gentamycin and attached to 35 mm dishes or multi-well plates (LifeSystemDesign) pre-coated with a laminin-5-rich extracellular matrix (28). For hormone secretion assays, approximately 15,000 cells were plated per dish, in three separated drops of 50 µl each. For bioluminescence recordings either 250 islets or approximately 50,000 separated cells were plated per well.

#### Quantitative RT-PCR (qRT-PCR)

Total RNA was prepared from homogenized islet cells using RNeasy® Plus Micro Kit (Qiagen). Ten nanograms of total RNA were reverse transcribed (PrimeScript RT reagent kit; Takara) and pre-amplified (TaqMan PreAmp Master Mix; Applied Biosystems) following the manufacturer's instructions. Specific target gene mRNA levels were analyzed by real-time quantitative PCR using the LightCycler technology (LC480; Roche Diagnostics). Mean values of gene expression were calculated from technical duplicates of each qRT-PCR analysis and normalized to the housekeeping gene *Hypoxanthine guanine phosphoribosyl transferase*  (*Hprt*) exhibiting no significant variability of its expression level throughout each experiment, and therefore served as internal control. Primers used for this study are listed in **Table 1**.

#### Hormone Secretion Measurements

Insulin and glucagon basal secretion assays were performed on approximately 15,000 attached separated α- or β-cells at 3 days *in vitro*. For basal insulin secretion assessment, cells were washed in KRB solution (Krebs-Ringer bicarbonate, pH 7.4, supplemented with 0.3% BSA) containing 2.8 mM glucose for 1 h, with subsequent incubation in KRB solution, containing 2.8 mM glucose for 30 min at 37°C in a cell culture incubator. For glucagon assessment, cells were washed in KRB solution containing 7 mM glucose for 1 h, followed by additional 30 min incubation in the same solution at 37°C in a cell culture incubator. To assess the islet cell acute secretory response, FACS-separated α- (*Venus*-positive) and β- (*Cherry*-positive) cells were subjected to 2 h incubation in KRB solution containing 5.6 mM glucose, followed by additional 30 min incubation in KRB containing 5.6 mM glucose (basal condition), 30 min incubation in KRB containing 5.6 mM glucose supplemented with 10 mM arginine (stimulated secretion for *Venus*-positive cells) or in KRB containing 16.7 mM glucose (stimulated secretion for *Cherry*-positive cells), and additional 30 min incubation in KRB containing 5.6 mM glucose (re-basal condition). At the end of each experiment, cells were lysed in Acid/Ethanol mixture (1.5% HCl/75% Ethanol) for glucagon or insulin residual contents measurements. Insulin and glucagon concentrations were quantified in the supernatants and in lysed cells by the Mouse Insulin or Glucagon ELISA kits (Mercodia). For acute secretion assays, glucagon and insulin values were expressed as a percentage of appropriate released hormone to residual cell content.

#### Islet Cells Synchronization and Circadian Bioluminescence Monitoring

Adherent islets/islet cells were synchronized by a 1-h pulse of forskolin 10 µM (Sigma) prior to continuous bioluminescence recording in RPMI, supplemented with 100 µM luciferin (NanoLight Technology) (20). Photon counts of each well were

integrated during 1 min, over 24 min intervals. For detrended time series, raw luminescence signals were smoothened by a moving average with a window of 24 h, allowing for a less biased comparison of bioluminescence values across experiments with regard to the measured circadian parameters (20).

#### RESULTS

### Separating Primary Mouse **α**- and **β**-Cells

In order to simultaneously label α- and β-cells within the islets, a *ProGcg*-*Venus*/RIP-*Cherry*/*Per2*:*Luc* reporter mouse line was established (**Figure 1**), allowing for the highly specific separation of nearly pure endocrine cell populations. To this end, following islet isolation and gentle trypsinization, dispersed cells were sorted by FACS based on cellular fluorescence characteristics, cell size (based on the data of Forward Scatter detector, FSC), and granularity (based on the data of Side Scatter detector, SSC; see **Figures 2A,B**). Of note, *Cherry*-positive cells showed higher cell granularity than *Venus*-positive cells (compare two histograms in **Figure 2B**). In parallel, cell viability was assessed by utilizing DRAQ7™, a dye that binds to DNA when cell membrane permeability is altered after initiation of cell death. Overall viability across preparations was near 90% (**Figure 2C**), with almost a 10-fold higher percentage of cell death for *Cherry*-positive cells (up to 20%) than for *Venus*-positive cells (**Figures 2D,E**), suggesting a higher sensitivity of *Cherry*-positive cells to the islet isolation, trypsinization, and/or sorting processes. The average number of harvested viable *Cherry*-positive cells per mouse was 39,292 ± 6,887, which is more than threefold higher than for *Venus*-positive cells (12,521 ± 1,885) (**Figure 2F**), reflecting the physiological ratio between these two cell types within mouse islets (22). In addition, the purity of the obtained *Venus*- and *Cherry*-positive cell populations was assessed by a complementary round of FACS analysis (**Figures 3A–G**). According to this analysis, the α-cell population comprises more than 90% viable *Venus*-positive cells without detectable contamination with *Cherry*-positive cells (**Figures 3A–C,G**), while the β-cell population contained up to 96% viable *Cherry*-positive cells without visible *Venus*-positive contaminants (**Figures 3D,E,G**). Finally, the morphological examination of sorted cell populations with a fluorescent microscope confirmed their high purity (**Figure 3H**).

### Primary Culture of Mouse **α**- and **β**-Cells

Insulin and glucagon transcript expression levels were assessed in separated *Venus*- and *Cherry*-positive cells by qRT-PCR analysis. As expected, *Gcg* transcription was the highest in the *Venus*-positive population, further confirming the α-cell identity, while both insulin transcripts *Ins1* and *Ins2* were abundant in *Cherry*-positive cells, indicating their β-cell identity (**Figure 4A**). Importantly, expression levels of the opposite cell hormone genes (*Ins1* and *Ins2* in *Venus*-positive cells, and *Gcg* in *Cherry*-positive cells) were more than 10-fold (for *Ins1* and *Ins2*) and 1,000-fold (for *Gcg*) lower, further confirming the satisfactory purity of the α- and β-cell populations. Furthermore, *Cherry*-positive cell population expressed β-cell-specific transcription factor *MafA* at the levels which were 1,000-fold higher compared to this transcript expression in *Venus*-positive counterparts (**Figure 4B**).

For cell culture, separated α- and β-cells were plated on plastic dishes covered with laminin-enriched matrix, which improves islet cell survival and function, and keeping them from de-diffentiation (28). Attached islet cells were maintained *in vitro* for at least 10 days. During the first 48 h of culture, attached β-cells formed monolayer cell aggregates resembling pseudo islets, while α- cells formed small domed structures composed by a few cells or stayed separated (**Figure 4C**). Immunofluorescence analysis demonstrated that *Venus*-positive cells in culture co-localized with glucagon-specific antibody, whereas *Cherry*-positive cells co-localized with insulinspecific antibody (**Figure 4D**), further validating these cell identity and high purity of α- and β-cell fractions.

Hormone secretion assays, performed after 3 days in culture, detected high basal levels of insulin in the supernatant of β-cells, and high basal levels of glucagon in α-cell supernatants, released during 30 min in the presence of constant glucose concentrations (2.8 mM for insulin and 7 mM for glucagon; **Figure 4E**). Noteworthy, secretion of the opposite hormones (glucagon by β-cells and insulin by α-cells) was below the detection level. Importantly, incubation of cultured α-cells with arginine induced about 2.5-fold increase in glucagon secretion, whereas incubation of cultured β-cells with high glucose stimulated secretion of insulin above two-fold (**Figure 4F**). Taken together, these data suggest that the here-described methodology ensures highly specific and efficient separation of the two main populations of islet cells, resulting in nearly pure populations of viable and functional α- and β-cells, which can be maintained in culture.

### Assessment of Circadian Oscillator Properties in Separated **α**- and **β**-Cells

In view of the complexity of the separation procedure by FACS, we next explored if separated α- and β-cell populations maintain their circadian properties, as was previously reported to be the case in isolated intact islets (11, 18). To this end, mRNA levels for selected core clock transcripts have been assessed in sorted islet cells isolated every 4 h during 24 h. The qRT-PCR analysis revealed pronounced rhythmic patterns for *Per1* and *Per2* genes over 24 h, while the oscillation of *Clock* transcription was shallow, in agreement with previous reports (**Figure 5**) (11, 19).

Moreover, we explored cell-autonomous molecular clocks in separated α- and β-cell primary cultures synchronized *in vitro*. Similar to intact islets (**Figure 6A**), populations of pure α-cells (**Figure 6B**) and β-cells (**Figure 6C**) responded to a 1 h synchronizing pulse of forskolin by demonstrating high-amplitude self-sustained circadian oscillations of *Per2:Luc* reporter expression for at least 5 consecutive days following synchronization. Synchronizing effect of forskolin on cultured α- and β-cells was specific, since medium change alone had little synchronizing effect on the *Per2:Luc* expression in both cell types (**Figures 6B,C**). Collectively, these data suggest that circadian oscillations persist in pure populations of α- and β-cells following islet isolation, trypsinization, and FACS separation procedures.

#### DISCUSSION

The major obstacle for studying α- and β-cells is that they are organized in the tight three-dimensional structure within the

FIGURE 2 | Fluorescence-activated cell sorting of pancreatic islet cells from ProGcg-Venus/RIP-Cherry/Per2:Luc mice. (A) Dot plots indicating good separation of *Venus*- and *Cherry*-positive cells from gated cell population. (B) Histograms representing number of analyzed events versus cell complexity (granularity) of the sorted cells. Note that *Cherry*-positive β-cells have greater granularity when compared to *Venus*-positive α-cells. (C) Dot plots for DRAQ7-based assessment of sorted cell viability, indicating approximately 90% alive cells in the preparation. (D) Representative histograms showing number of dead α- and β-cells among analyzed events and (E) corresponding quantitative data from 12 independent experiments (three to six mice per experiment), suggesting higher cell death for β-cell population, as compared to α-cells (paired, two-tailed Student's *t*-test). (F) Histogram indicating average numbers of obtained α- and β-cells per mouse (*N* = 12 independent experiments with three to six mice per experiment), statistical difference illustrate more than three times greater content of β-cells within sorted population, as compared to α-cells (paired, two-tailed Student's *t*-test). Data expressed as mean ± SEM.

FIGURE 3 | Purity of **α**- and **β**-cell populations, separated by fluorescence-activated cell sorting (FACS) form ProGcg-Venus/RIP-Cherry/Per2:Luc mouse islets. Sorted samples of *Venus*- and *Cherry*-positive cells were subsequently analyzed by the second run of FACS. (A–C) FACS-assessed purity analysis of obtained α-cell population: (A) dot plot presenting distribution of cell size [forward scatter detector (FSC), *x*-axis] versus granularity [side scatter detector (SSC), *y*-axis] of *Venus*-positive cells within gated cell population; (B) dot plot showing sample enrichment for *Venus*-positive cells and complete absence of *Cherry*-positive cells; (C) DAPI-assessed analysis of α-cell viability after sorting, indicating 97% viable cells in the sample. (D–F) FACS-assessed purity analysis of obtained β-cell population: (D) dot plot presenting distribution of cell size (FSC, *x*-axis) versus granularity (SSC, *y*-axis) of *Cherry*-positive cells within gated cell population; (E) dot plot showing enrichment of sample for *Cherry*-positive cells and absence of *Venus*-positive contaminants; (F) DAPI-assessed analysis of β-cell viability after sorting, indicating approximately 97% viable cells in the sample. (G) Histogram representing average purity in analyzed samples (data are presented as mean ± SEM, *N* = 6 experiments for α-cells; and *N* = 10 for β-cells). (H) Representative fluorescent images illustrating sorted *Venus*- and *Cherry*-positive cell populations (α- and β-cells, respectively). Scale bar = 100 µm.

pancreatic islet. We successfully overcame this problem by utilizing transgenic mice, specifically expressing the *ProGcg*-*Venus*  reporter in α-cells (23) and the RIP-*Cherry* reporter in β-cells (24) (**Figure 1**), allowing to separate these two cell populations by FACS with high viability and purity (**Figures 2** and **3**). The here-described methodology allows for extracting up to 40,000 β-cells and 12,500 α-cells per mouse (**Figure 2F**), and culturing thus separated primary α- and β-cells for at least 10 days. Importantly, the cell ratio after sorting reflects proportion between α- and β-cells in mouse pancreatic islets *in vivo* (60–80 versus 15–20%, respectively) (22, 29). The value obtained by FACS for this islet cell ratio following isolation and separation provides an estimation for the islet cell composition, which gives an advantage for studying the altered ratio between α- and β-cells upon different pathological conditions like obesity, T2D and others.

#### FIGURE 4 | Continued

Endocrine identity of sorted *Venus*- and *Cherry*-positive cell populations, isolated from ProGcg-Venus/RIP-Cherry/Per2:Luc mouse islets by fluorescence-activated cell sorting (FACS). (A,B) Quantitative RT-PCR analysis of (A) insulin (*Ins1* and *Isn2*) and glucagon (*Gcg*) transcripts; or (B) V-Maf avian musculoaponeurotic fibrosarcoma oncogene homolog A (*MafA*) β-cell-specific transcript in freshly isolated *Venus*- and *Cherry*-positive cells. Data are presented as relative expression of the target gene normalized to the *Hprt* housekeeping gene (mean ± SEM, *N* = 3 independent experiments for (A) and *N* = 4 experiments for (B) with six mice each). (C) Representative bright-field images of cultured *Venus-* and *Cherry*-positive cells 48 h after plating. (D) Representative fluorescent images of ProGcg-Venus and Insulin-Cherry transgenic proteins (left panels) colocalizing with glucagon and insulin immunostaining, respectively (merge, right panels). Analysis of purified *Venus-* and *Cherry*-positive cell cultures prior the immunostaining showed no cross-contaminations. No *Venus*-positive cells appeared positive for insulin, and no *Cherry*-positive cells appeared positive for glucagon by immunostaining. Scale bars = 100 µm. (E) FACS-separated *Venus*- and *Cherry*-positive cells were plated in 35 mm dishes (approximately 15,000 cells per dish, as described in Section "Material and Methods"), and subjected to 30 min incubation with 2.8 mM glucose (for basal insulin assessment) or with 7 mM glucose (for basal glucagon assessment). Hormone levels were assessed in the incubation solution with Mouse Insulin and Glucagon ELISA kits (Mercodia). Data are expressed as absolute values (mean ± SD) for *N* = 6 mice used for two independent islet preparations. (F) Acute secretion assays performed on cultured *Venus*- and *Cherry*-positive cells 48 h after plating. Cell supernatants were collected for three successive 30-min periods (30 min basal condition, 30 min stimulated condition, and 30 min re-basal condition, as described in Section "Material and Methods") after 2 h of cell depletion in KRB solution. Hormone levels were assessed with Mouse Insulin and Glucagon ELISA kits and normalized to the residual hormone content in the end of the experiment. Data are expressed as percentage of glucagon or insulin from cell contents (mean ± SD; *N* = 2 experiments for *Venus-*positive cells and *N* = 4 experiments for *Cherry*-positive cells (two mice per experiment); \* *p* < 0.05, two-tailed Student's *t*-test).

corresponds to 07:00 a.m., and ZT12 corresponds to 07:00 p.m.), with constant access to the drinking water and night-restricted feeding regime for 2 weeks prior to islet collection. Pancreatic islets were isolated around-the-clock every 4 h and subjected to the FACS separation procedure. Temporal profiles of selected core clock transcripts (*Clock*, *Per1*, and *Per2*) were assessed in freshly extracted mRNA at each time point from α- and β-cell populations by quantitative RT-PCR analysis. Transcript levels were normalized to *Hprt* expression. Data expressed as mean ± SD (*N* = 2 biological replicates of 6 mice at each time point).