Gene Expression Profiling in the Hibernating Primate, Cheirogaleus Medius.

Hibernation is a complex physiological response that some mammalian species employ to evade energetic demands. Previous work in mammalian hibernators suggests that hibernation is activated not by a set of genes unique to hibernators, but by differential expression of genes that are present in all mammals. This question of universal genetic mechanisms requires further investigation and can only be tested through additional investigations of phylogenetically dispersed species. To explore this question, we use RNA-Seq to investigate gene expression dynamics as they relate to the varying physiological states experienced throughout the year in a group of primate hibernators-Madagascar's dwarf lemurs (genus Cheirogaleus). In a novel experimental approach, we use longitudinal sampling of biological tissues as a method for capturing gene expression profiles from the same individuals throughout their annual hibernation cycle. We identify 90 candidate genes that have variable expression patterns when comparing two active states (Active 1 and Active 2) with a torpor state. These include genes that are involved in metabolic pathways, feeding behavior, and circadian rhythms, as might be expected to correlate with seasonal physiological state changes. The identified genes appear to be critical for maintaining the health of an animal that undergoes prolonged periods of metabolic depression concurrent with the hibernation phenotype. By focusing on these differentially expressed genes in dwarf lemurs, we compare gene expression patterns in previously studied mammalian hibernators. Additionally, by employing evolutionary rate analysis, we find that hibernation-related genes do not evolve under positive selection in hibernating species relative to nonhibernators.

Introduction

To survive in highly seasonal environments, mammals employ a range of phenotypic and behavioral tactics to cope with fluctuating energetic demands. Hibernation is among the most extreme of these strategies. This seasonally expressed period of heterothermy represents a radical deviation from mammalian homeostasis, involving complex physiological and behavioral changes that conserve up to 90% of the energy required to remain euthermic during the winter (Geiser and Heldmaier 1995; Morin and Storey 2009). Classical hibernation behavior is defined by periods of torpor (i.e., controlled reductions in metabolic rate, body temperature [Tb], and heart rate) and metabolically active periods of rewarming known as interbout arousals (IBAs).

The remarkably widespread phylogenetic diversity of hibernation among mammalian lineages raises the question: “Are there fundamental similarities in the molecular architecture underlying the hibernation phenotype?” (Srere et al. 1992; Carey et al. 2003; Villanueva-Cañas et al. 2014). Hibernators are found in 14 orders of mammals that range from arctic, temperate, and tropical regions, and represent all three of the deepest branches of the Class Mammalia: Monotremata, Marsupialia, and Placentalia (Srere et al. 1992; Carey et al. 2003; Andrews 2007; fig. 1). The current understanding is that the hibernation phenotype arises due to selective expression of genes that are correlated with orchestrating physiological changes that allow a hibernator to experience recurrent torpor bouts during the winter (Srere et al. 1992; Wilson et al. 1992; Mominoki 1998; Bauer et al. 2001; Buck et al. 2002; Epperson and Martin 2002; Eddy 2004; Mominoki et al. 2005; Williams 2005; Eddy et al. 2006; Yan 2006; Suozzi et al. 2009; Epperson et al. 2010; Shao et al. 2010; Fedorov et al. 2011; Hampton et al. 2011; Chow et al. 2013; Hampton et al. 2013; Schwartz et al. 2013; Seim et al. 2013; Emirbekov and Pashaeva 2014; Lei et al. 2014; Grabek et al. 2015; Vermillion et al. 2015). However, it is unknown whether distantly related species use the same key pathways for activating and maintaining the hibernation phenotype. Our study contributes to a broader sampling across the mammalian phylogeny by incorporating a group of primate hibernators, the dwarf lemurs of Madagascar (genus Cheirogaleus), and thus has the potential to provide unique insight into the question of a universal genetic regulatory pathway of hibernation.

. 1.—

The patchwork distribution of heterothermic species within the mammalian phylogeny demonstrates that hibernation is likely a retained ancestral trait. Orange box denotes Monotremata, yellow box indicates Marsupialia, and green box defines Placentalia. Families colored in blue text indicate the inclusion of at least one species that is heterothermic. Representative heterothermic species are pictured on right. No distinction was made between summer heterothermy and winter heterothermy. Phylogeny modified from dos Reis et al. (2012).

Dwarf lemurs belong to one of two clades within the primate lineage to manifest depressed metabolic rates through torpor. The closely related mouse lemurs, which, such as dwarf lemurs, are members of the family Cheirogaleidae, are capable of torpor (Ortmann et al. 1997; Schmid 2001), although torpor bouts in these animals are typically shorter than 24 h. The other primate clade to express torpor is the Lorisiformes. The Asian pygmy slow loris was most recently documented to demonstrate bouts of sustained metabolic suppression under semicaptive conditions (Ruf et al. 2015), whereas anecdotal reports of torpor have also been reported for their sister lineage, the bushbabies, under extreme energetic demands (e.g., Nowack et al. 2013). In the western regions of Madagascar, fat-tailed dwarf lemurs (Cheirogaleus medius) engage in hibernation for up to 7 months wherein body temperature approximates ambient temperature (Dausmann et al. 2004, 2005, 2009). Madagascar is highly seasonal, with severe shortages of rain and associated resources during the winter. In response, fat-tailed dwarf lemurs enter hibernation to reduce energetic costs when dry spells occur (Dausmann et al. 2005, 2004, 2009). When occupying well-insulated hibernacula, fat-tailed dwarf lemurs display a classic hibernation pattern including characteristic IBAs; however, when hibernating in poorly insulated tree holes, they can abandon an active defense of Tb and switch to a poikilothermic state where Tb thermoconforms to ambient without the demand for IBAs (Dausmann et al. 2004). This is in contrast to every other known hibernating species (except members of the family Ursidae and Tenrecidae; Toien et al. 2011, Lovegrove et al. 2014). Presumably the ability to thermoconform under poorly insulated conditions allows these tropical primates to engage in hibernation, even at widely fluctuating (i.e., oscillations of more than 30° in a single day), elevated Ta (Dausmann et al. 2009, 2005).

Dwarf lemurs can more than double their body weight by building up reserves of white adipose tissue (WAT) preceding hibernation and then use these stores as an energy reservoir during the prolonged period of physical inertia (Fietz and Dausmann 2007). In captivity, even in the absence of environmental triggers such as food scarcity and lowered temperatures, fat-tailed dwarf lemurs nonetheless exhibit seasonal body mass cycles. This seasonal response, irrespective of environmental conditions, strongly suggests that underlying clock-like molecular mechanisms are responsible for seasonal fat deposition (Perret et al. 1998; Perret and Aujard 2001; Dark 2005).

WAT is an especially informative organ for study in the context of the hibernation phenotype given that fat deposition and metabolism are critical aspects of survival during hibernation. Fat-storing hibernators switch their metabolic economy from carbohydrates during the active season, to lipids during the hibernation season as a means to conserve energy (Carey et al. 2003; Boyer and Barnes 1999; Buck et al. 2002). Investigations of how annual body mass cycles are governed by underlying molecular mechanisms in WAT have been conducted in multiple mammalian lineages using candidate gene approaches (Wilson et al. 1992; Herminghuysen et al. 1995; Boyer et al. 1998; Bauer et al. 2001; Buck et al. 2002; Demas et al. 2002; Eddy 2004; Kabine et al. 2004; Eddy et al. 2005), and recently using next-generation sequencing (Hampton et al. 2011) in a subset of mammalian model hibernators. These studies have laid the groundwork for establishing the paradigm that differential gene expression is correlated with physiological changes throughout the year. Though one previous study employed a longitudinal approach to investigate protein levels in the blood serum of black bears (Chow et al. 2013), ours is the first to examine levels of gene expression in biological tissues (i.e., WAT) from the same animals at multiple times points during the year. This is a powerful approach for controlling for interindividual variation, which might otherwise confound differential gene expression analyses.

It is currently unknown how seasonal body mass cycles are regulated in dwarf lemurs, and if dwarf lemurs conform to the conventional paradigm that changes in metabolic economy occur in conjunction with variations in gene expression throughout the course of the circannual cycle. Therefore, in the present study, we use gene expression data generated from an RNA-Seq experiment, sampled longitudinally in multiple animals, to examine the genetic controls behind the metabolic switch from carbohydrate to lipid metabolism in dwarf lemurs. We compare our results with previously published data from model mammalian species to provide further insight into the question of whether the hibernation phenotype in distantly related mammalian hibernators is controlled by the same regulatory pathways (Villanueva-Cañas et al. 2014).

Materials and Methods

Animals and Care Regimen

All animals examined in this study were captive born and have been maintained under captive conditions throughout their lifespan. The fat-tailed dwarf lemur colony maintained at the Duke Lemur Center (DLC) offers a unique resource for systematically sampling key biological tissues under carefully controlled conditions using a longitudinal sampling approach. A total of four (two male and two female) fat-tailed dwarf lemurs were used throughout this study with ages ranging from 8 to 13 years old. At the beginning of the hibernation season at DLC (November), dwarf lemurs were moved into a specially constructed hibernaculum, which allowed independent control of photoperiod and temperature. Animals were caged separately, but allowed visual and auditory access to other individuals. Individuals were provided water ad libitum and a standard DLC diet of Purina primate chow, fruits/vegetables, and mealworms. Under captive conditions, dwarf lemurs are provided multiple nesting tubes to mimic their natural habitat, where they hibernate in tree holes, and are therefore still exposed to photoperiodic triggers. Study animals are exposed to simulated natural conditions, including a 66% reduction in caloric input preceding torpor to simulate reduced food resources, shortened daylight hours, and lowered ambient temperature, so as to induce seasonal torpor during the months of January and February.

Experimental Collection Points

Three experimental time points were used: October (Active 1), December (Active 2), and February (Torpor) to represent a wide range of physiological changes that coincide with changes in seasonality (fig. 2). Physiological state at each collection point was verified by rectal Tb and body mass was recorded at the time of sampling (table 1; individual animal’s data can be found in supplementary table S1, Supplementary Material online). The only difference between October and December was a gradual decrease in daylight hours from 11.5 to 9.5 h of daylight, simulating the transition from summer into fall. Our initial goal was to compare three distinct physiological states, including a period of extreme fattening during December (Active 2). Due to federally mandated standards, however, animals included in our study were not allowed to feed ad libitum during the fall months which, in wild populations, is crucial for maximum fat accumulation preceding the hibernation season. This dietary restriction may therefore have had important consequences for the comparison of Active 1 relative to the Active 2 collection points.

. 2.—

Schematic representation of timeline for experimental design. Solid black line indicates cyclic body weight changes during the year. Dotted red line depicts ambient temperatures used during study period. Circles along the x axis indicate photoperiodic changes that occurred during the study period, mimicking the switch from autumn to winter day lengths.

Table 1

Collection Points and Corresponding States of Activity, Mean Body Mass (±SD), Mean Body Temperature (Rectal Tb; ±SD), Ambient Temperature (Ta), and Photoperiod during Sample Collection

Collection Point	Physiological State	Body Mass (g)	Tb (°C)	Ta (°C)	Photoperiod
October 2012	Active 1	238.7 ± 32.1	33.9 ± 0.1	25	LD 11.5:12.5
December 2012	Active 2	237.3 ± 29.3	35.3 ± 1.8	25	LD 9.5:14.5
February 2013	Torpor	196.0 ± 57.2	19.8 ± 4.1	15	LD 10.5:13.5

During a torpor bout, dwarf lemurs exhibit decreases in Tb to as low as 15 °C, reductions in heart rate to 3–5 beats per minute, and episodes of irregular and infrequent breathing (Krystal et al. 2013; Dausmann et al. 2004). We sampled WAT from animals when they were 2–3 days into a torpor bout as determined by an absence of feeding.

Tissue Collection

We successfully used minimally invasive sampling techniques to acquire tissue samples of WAT. Animals were anesthetized using Ketamine (10 mg/kg; Ketaject, Bioniche Teoranta, Galway, Ireland) and Midazolam (0.25 mg/kg; Bedford Laboratories, Bedford, OH) by a DLC veterinarian. Sampling procedures took place at a consistent daily time point (between 9:00 and 10:00 AM) to account for any circadian fluctuations influencing gene expression profiles. Tissue samples of WAT were obtained using needle biopsy from the section of the tail where fat storage was most concentrated (Fietz et al. 2003). The area to be sampled was shaved and skin was cleaned with chlorhexidine solution followed by alcohol three alternating times. Skin was held taught and a 16-gauge tru-cut biopsy needle (Jorgenson Laboratories, Loveland, CO) was inserted in target area. Biopsy was removed and tissue sample was expelled directly into RNA stabilization solution (Qiagen, Valencia, CA). Samples were placed in 4 °C overnight, and the next morning moved to storage in −20 °C until further processing. The 4 °C incubation prior to freezing was done to allow cross-comparisons with samples collected from wild individuals for a future investigation.

RNA Purification and Whole Transcriptome Amplification

The entire tissue sample (15–20 mg) was used for RNA isolation procedure. All RNA extractions were conducted on the same day after tissue sampling was complete by the lead author (S.L.F.) to minimize intersample variability. Total RNA was purified from WAT using an optimized TRIzol protocol in conjunction with the Microarray Tissue Kit (Qiagen). Briefly, frozen tissues were homogenized in 500 μl TRI Reagent (Ambion, Grand Island, NY) using a hand-held rotor-stator homogenizer to provide efficient disruption and homogenization of samples. BCP (1-bromo-3-chloropropane) extraction was performed using 50 μl of BCP, followed by centrifugation at 12,000 × g for 15 min at 4 °C. Aqueous layer containing total RNA was transferred to a fresh 1.5-ml centrifuge tube and one volume of 70% EtOH was added. The entire sample volume was loaded onto Qiagen RNeasy filter columns and kit protocol was followed according to manufacturer’s instructions, including an on-column DNase step to remove any residual contaminating DNA. RNA integrity and concentration was assessed using the Agilent 2100 Bioanalyzer (Santa Clara, CA) at Duke University’s Institute for Genome Science and Policy’s Microarray Facility. RNA Integrity Number (RIN) values for our total RNA extractions ranged from 3.5 to 6.7 and total RNA concentrations ranged from 2 to 24 ng/μl. WAT has a high lipid content to nuclear content ratio and extracting sufficient total RNA for downstream RNA sequencing is therefore problematic. To that end, we completed a whole transcriptome amplification step on all total RNA extractions prior to Illumina sequencing.

To account for the wide range in extracted RNA concentrations, purified RNA was diluted to 2 ng/μl in a 5 μl total volume and then subjected to whole transcriptome amplification using NuGEN’s Ovation RNA-Seq V2 kit (San Carlos, CA) according to manufacturer’s instructions. This kit provides a rapid method for preparing amplified cDNA from total RNA for downstream RNA-Seq applications. It employs a single primer isothermal amplification method to amplify total RNA into double-stranded cDNA and depletes rRNA without preselecting mRNA. Amplified cDNA samples were then purified using the MinElute Reaction Cleanup Kit (Qiagen) according to manufacturer’s protocol. cDNA concentrations postamplification ranged from 347 to 461 ng/μl. In a preliminary analysis using unamplified versus amplified RNA extracted from rat WAT, we have confirmed that whole transcriptome amplification does not introduce significant bias in relative mRNA frequencies (Faherty et al. 2015). Amplified cDNA samples from dwarf lemurs were then sent to Duke University’s Genome Sequencing Shared Resource for library preparation and sequencing.

Library Preparation and Illumina Sequencing

Previously amplified cDNA libraries were prepared for sequencing using Illumina’s TruSeq DNA Sample Preparation Kit. Final library size distribution was determined using Agilent Bioanalyzer 2100 and insert size ranged from 206 to 246 bp. Libraries were pooled and sequenced on two lanes of the Illumina HiSeq2000 platform (San Diego, CA) using the rapid-run mode with 150-bp paired-end reads. The content of each library was divided by half and each half sequenced on one of the two lanes. This was done to avoid lane-related batch effects. A summary of the RNA-Seq data is displayed in supplementary table S2, Supplementary Material online.

Quality Control and Filtering of Sequencing Reads

A schematic of the computational pipeline for assembling C. medius transcripts from raw sequencing reads, calculating expression levels, and determining differentially expressed (DE) genes between conditions is shown in figure 3.

. 3.—

Overview of the bioinformatics pipeline used in this study.

Quality control processing of the RNA-Seq data was completed using FastQC software developed by Babraham Bioinformatics (Andrews 2010). This software calculates global statistics based on Phred scores and nucleotide composition, including plots of positional Phred scores across the length of all reads, positional nucleotide composition, and relative k-mer enrichment over read length which was used to define parameters for downstream filtering using the software Trimmomatic (Bolger et al. 2014). We removed the Illumina adaptor sequences and trimmed the ends of the reads by checking the first and last three bases and retaining reads with a Phred score >30, with a required minimum length of 100 bp. The median trimmed length for the forward strand was 151 ± 13.34 nt, whereas the reverse strand had a median trimmed length of 151 ± 9.50 nt.

De Novo Transcriptome Assembly

We constructed a reference de novo assembly by pooling the RNA-Seq reads from all samples and using the bioinformatic pipeline Trinity (version: 20130814) (Haas et al. 2013) on a computing cluster with 48 cores and 512 Gb RAM. We used the following command line: Trinity.pl –seqType fq –JM 90G –left all_1.fastq –right all_2.fastq –output assembly_output –CPU 22. We obtained 1,295,606 putative transcripts corresponding to 776,170 gene models with an N50 of 981 nucleotides. The number of initially assembled transcripts was very high, due to inclusion of any assembled transcripts longer than 200 nucleotides. Short transcripts are typically discarded downstream. Subsequently, we filtered out putative artifacts, due to polymorphisms and/or sequencing errors, with the program CD-HIT-EST (Fu et al. 2012). By employing a 95% similarity threshold the number of transcripts decreased to 1,206,443. Although the majority of transcripts was quite short (supplementary fig. S1, Supplementary Material online), we were able to reconstruct a large number of long transcripts (225,744 transcripts > 1 kb). Next, we aligned all the reads to the assembly using Bowtie2 (Langmead and Salzberg 2012) and used the Integrated Genome Viewer (IGV) (Thorvaldsdottir et al. 2013) for data visualization. Finally, we performed sequence similarity searches of the reconstructed transcripts against the human proteome (Ensembl v.73) with BLASTX 2.2.28+ (Camacho et al. 2009). Using an e-value < 10−4 we found 158,276 transcripts in C. medius with a detectable homologous protein in human. After eliminating redundant hits we detected 9,540 transcripts with homology to human proteins (>70% of the protein length). For differential gene expression analysis, we used the 15,420 transcripts with homology in humans that expressed in at least four samples. Details on the transcript assembly and comparisons to other protocols that were applied to the same data set can be found in supplementary figures S1 and S2 and supplementary table S3, Supplementary Material online.

Investigating the Influence of RNA Quality on Transcriptome Assembly

We investigate how RIN scores, which measure RNA integrity, affected the number of assembled transcripts and coverage at the 5′- and 3′-transcript ends. As expected, the number of reconstructed transcripts was positively correlated with the RIN score of the sample (supplementary fig. S2A, Supplementary Material online). However, we also noted that even in samples with very low RIN scores (≤5), the number of transcripts recovered was quite large (on average, 67.08% of the number recovered in the best sample).

Using an in-house python script, we clustered together all the reads mapping to each transcript and classified them as 5′ if they were located in the first 10% of the transcript, and 3′ if they were located in the last 10% of the transcript. When comparing the degree of enrichment or depletion of sequencing reads in either the 5′- or 3′-regions, we detected a statistically significant impoverishment in the 3′-region in all the samples (one sample Wilcoxon test, µ = 0.1). The number of reads in the 3′-region was about half the number expected if the reads were homogeneously distributed. Surprisingly, however, we found no correlation between the level of 3′-read depletion and the RIN score, suggesting that this bias is intrinsic to the RNA-Seq technique and largely independent of sample quality as indicated by RIN score (supplementary fig.S2B, Supplementary Material online).

Identification of Human Genes Expressed in Adipose Tissue Using GTEx Data

We downloaded gene expression data for the two RNA-Seq human adipose tissue samples available from the GTEx Consortium (GTEX-N7MS-0326-SM-4E3K2 and GTEX-NFK9-0326-SM-3MJGV; Lonsdale et al. 2013). We determined that 11,875 protein-coding genes were expressed at Fragments Per Kilobase of transcript per Million (FPKM) >0.5 in both samples, 8,682 of which corresponded to dwarf lemur assembled genes.[AQ6]

Estimating Transcript Abundance and Differential Gene Expression

To estimate transcript abundance, we used the software package RSEM (Li and Dewey 2011). We generated a matrix of normalized counts comprised the number of reads for each sample (columns) mapping to each reconstructed transcript (rows). To identify DE transcripts, we used the edgeR Bioconductor package (Robinson et al. 2009). edgeR is based on an overdispersed Poisson model that accounts for technical and biological variability. We used the following general linear model: ∼condition + RIN, to avoid any putative biases created by the range of RIN scores in our samples. We only tested C. medius genes with a significant hit against the human proteome (Ensembl v.73) using an e-value threshold of 10−4. We discarded any genes in which at least four samples had less than 0.5 counts per million. This was done to remove the effect of poorly supported transcripts in some samples. We identified DE genes for each sample collection point by analyzing data from three individuals (Tanager, Quetzal, and Osprey; see supplementary table S2, Supplementary Material online). One individual was removed (Towhee) because she lost a significant amount of weight prior to entering torpor, probably as a consequence of the food-restriction guidelines described in the section “Experimental Collection Points.” Our analyses were performed using pairwise comparisons, as we were focused on detecting differences between physiological states. The individual samples were used to assess the biological variability within each condition.

We corrected P values for multiple comparisons using the Benjamini–Hochberg method (Benjamini and Hochberg 1995) providing a P value cut-off for significance which controlled the false discovery rate (FDR) at 0.05. We performed a functional annotation clustering analysis of the DE genes using the Panther Classification System (Panther v.10) (Huaiyu et al. 2016). We interrogated groups of genes that were up- or downregulated during the Torpor time point in relation to the Active state.

Protein Evolutionary Rates in Hibernating Species

We extracted the protein sequence from the Trinity reconstructed transcripts in dwarf lemurs with the best BLASTX match against the human proteome. We selected cases in which we could detect one-to-one orthologous genes in ten species from Ensembl (v.82), three of which are capable of hibernation (marked with an asterisk): Homo sapiens, Bos taurus, Cavia porcellus, Mus musculus, Myotis lucifugus*, Pteropus vampyrus, Oryctolagus cuniculus, Otolemur garnettii, Spermophilus tridecemlineatus*, and Cheirogaleus medius*.

We generated multiple sequence alignments with the software Prank + F using the algorithm PALO for the selection of the most appropriate protein for each gene, instead of the longest one. This methodology has been shown to reduce false positives in the estimation of the rate of nonsynonymous to synonymous substitutions (dN/dS = ω) using CodeML (Fletcher and Yang 2010; Villanueva-Cañas et al. 2013). The alignments were manually trimmed afterwards. We compared two different models using CodeML in the PAML package (Yang 2007). The tree with branch distances for CodeML input was modified from Meredith et al. (2011). In the first model, leaves of hibernating species and nonhibernating species were allowed different ω. In the second model, all leaves had the same ω value. The likelihood ratio between the two models was obtained with LR = 2*(lnL1 − lnL2) with 1 degree of freedom.

Results

We performed de novo transcript assembly using Trinity by concatenating the reads from all 12 collected samples (supplementary tables S1 and S2, Supplementary Material online). Sequence similarity searches with BLASTX revealed that 158,276 transcripts had a detectable homolog in humans. Of these, some of the transcripts were partially reconstructed and were mapped to a redundant human protein. After clustering redundant hits, we annotated 9,540 unique C. medius genes with homology to human proteins. A substantial majority of them (8,682; ∼91%) were expressed in the two available human GTEx samples from adipose tissue.

The preservation of RNA in biological samples is especially challenging in studies performed in the wild, or when samples are problematic to acquire, such as biopsied tissue (Gallego Romero et al. 2014). To investigate the effect of the RIN scores on transcript assembly, we compared the number of detected transcripts in each sample with the range of RIN values. Although there was some effect across samples, even for samples with very low RIN values, we reconstructed a high number of transcripts with human homologs (supplementary fig. S2A, Supplementary Material online). The number of reads falling into 5′- and 3′-ends of transcripts appeared to be independent of the RIN value (supplementary fig. S2B, Supplementary Material online). This suggests that the effects are relatively small. Nevertheless, given that other studies have suggested that this kind of analysis is sensitive to RNA integrity (e.g., Gallego Romero et al. 2014), we explicitly took the RIN values into account when performing the differential gene expression analysis.

Identification of DE Genes

We estimated gene expression levels by mapping back the reads of each individual sample to the dwarf lemur transcript. Subsequently, we performed differential gene expression analysis between pairs of conditions with edgeR. Despite our sampling being longitudinal in nature, we were only interested in comparing each time point in a pairwise fashion. The linear model employed included the RIN value of the sample. Table 2 shows a summary of the number of DE genes (FDR < 0.05) uncovered. Only one gene was differentially expressed in the Active 1 versus Active 2 comparison (table 2). In contrast, more than 100 genes were differentially expressed in the Active 1 versus Torpor, and Active 2 versus Torpor pairwise comparisons, showing an overlap of 71%. This corresponds to 90 DE genes that consistently change expression patterns during Torpor relative to Active states.

Table 2

Total Number of Genes Assessed for Differential Expression and DE Genes among Different Physiological Conditions

Pairwise Comparison	DE Genes with Sequence Similarity to Human (e-Value < 10−4)
Active 1 versus Torpor	283
Active 2 versus Torpor	126
Active 1 + Active 2 versus Torpor	90
Active 1 versus Active 2	1

DE Gene Functions and Pathways

Functional enrichment analysis was performed with the Panther Classification System after grouping genes both upregulated and downregulated during Torpor as distinct groups. This analysis did not show any clear enrichment in DE genes after correcting for multiple testing; however, given the manageable number of DE genes, we examined the function of each gene individually. We highlight genes that show the same pattern of differential expression across all individuals when comparing each physiological state. The complete list of DE genes for all time point comparisons is provided in supplementary file S1, Supplementary Material online. Based upon body weight and core Tb data, we do not discriminate between Active 1 and Active 2, but instead treat them as biological replicates to characterize gene expression profiles during an Active physiological state. When comparing Active to Torpor states, we find genes involved in metabolic processes, DNA and protein binding, cellular adhesion, and blood coagulation/circulation to be differentially expressed in a seasonal pattern.

Our experimental design also allowed for investigations into the influence of photoperiod on gene expression profiles. Between Active 1 and Active 2 time points, photoperiod was reduced to simulate the transition from summer to fall. We compared each Active collection point individually with our Torpor collection point. Genes involved in metabolic pathways and fat storage that displayed upregulation during Torpor include adipose acyl-CoA synthase-1 (ACSL1; 2.53-fold over Active 2), aminoadipate-semialdehyde synthase (AASS; 4.54-average fold over combined Active 1 and Active 2), pyruvate dehydrogenase kinase 4 (PDK4; 1.66-fold over Active 1), and phosphofructokinase platelet (PFKP; 3.29-average fold over combined Active 1 and Active 2) (fig. 4).

. 4.—

Boxplot graph showing significant expression (BH < 0.05) changes between physiological states for highlighted genes mentioned in the text. Y axis shows the number of reads in log10 scale that map to each reconstructed gene, whereas boxplot whiskers show the range of reads between individuals.

We also find other genes involved in feeding behavior to have seasonal expression patterns. We find leptin (LEP) and its receptor, LEPR, to be downregulated during Torpor when compared with the Active 1 state (decreased 3.72-fold and 2.31, respectively). Genes implicated in the mammalian circadian clock are also differentially expressed. During Torpor, two isoforms of nuclear factor, interleukin 3 regulated (NFIL3), a transcription factor that regulates Period clock genes (PER1 and PER2), are upregulated 2.57-fold and 2.21-fold over Active 1. PER1 is upregulated 1.85-fold over Active 1. Genes involved in the blood coagulation cascade show variable expression patterns throughout the year. When comparing Torpor to the Active 1 state we find that a member of the annexin family (ANXA4; 1.80-fold decrease during Torpor) and coagulation factor VIII (F8; 2.50-fold decrease) display differential gene expression (fig. 4).

Extreme Changes in Gene Expression across Collection Points

A number of genes were found to display notable seasonally dependent changes, as determined by the average log fold change between the combined Active 1 and Active 2 versus Torpor pairwise comparison. The most drastic examples include desmoplakin (DSP; 10.04-fold), membrane bound O-acyltransferase domain containing 2 (MBOAT2; 9.23-fold), exophilin 5 (EXPH5; 8.09-fold), desmocollin 3 (DSC3; 7.58-fold), haptoglobin (HP; 7.22-fold), and hydroxysteroid (17-beta) dehydrogenase 6 (HSD17B6; 6.29-fold)

Shared Genes among Mammalian Hibernators

Using gene expression profiles generated from our RNA-Seq experiment for the combined Active 1 and Active 2 versus Torpor pairwise comparison, we used previously published data sets from multiple ground squirrel species, the American black bear, and bats to look at common genes identified as differentially expressed in conjunction with changes in physiology. Table 3 highlights the DE genes we identify in dwarf lemurs that are also identified as differentially expressed in other hibernating species.

Table 3

Differentially Expressed Genes in Active 1 + Active 2 versus Torpor Collection Point with Homolog in Another Hibernating Species Determined from Previous Studies

Gene Name	Symbol	logFC	References
Metabolism and metabolic processes
Aminoadipate-semialdehyde synthase	AASS	5.02	Epperson et al. 2010; Shao et al. 2010
Phosphofructokinase, platelet	PFKP	3.66	Lei et al. 2014
Basic cellular processes
Inhibitor of DNA binding 4	ID4	−2.42	Schwartz et al. 2013
Phosphodiesterase 1A, calmodulin-dependent	PDE1A	−3.30	Schwartz et al. 2013
Pleckstrin homology-like domain, family B, member 2	PHLDB2	−1.48	Fedorov et al. 2011
S100 calcium binding protein A1	S100A1	−2.22	Schwartz et al. 2013; Emiberkov and Pashaeva 2014; Vermillion et al. 2015
Circulation and blood coagulation
Haptoglobin	HP	7.06	Mominoki 1998; Mominoki et al. 2005; Yan et al. 2008; Chow et al. 2013; Vermillion et al. 2015

Protein Evolutionary Rates

The assembly of a transcriptome de novo gave us the opportunity to perform protein sequence comparisons between the dwarf lemur and other species. In particular, we were interested in testing whether hibernation-associated genes evolved faster in hibernating species versus in nonhibernating species. We selected a set of 12 genes that had one-to-one orthologs in nine other mammals, including two other hibernating and seven nonhibernating species (see Materials and Methods for the list of species). We estimated the rate of nonsynonymous versus synonymous substitutions (dN/dS) and compared the likelihood of two models, one in which the dN/dS was independent of the hibernating status and another one in which the dN/dS was assumed to be different between hibernating and nonhibernating species (supplementary fig. S3, Supplementary Material online). In only 2 of 12 cases the second model was more plausible (P value < 10 − 3; supplementary file S2, Supplementary Material online). Although the gene showed a significant increase in dN/dS in the hibernation-associated branches, dN/dS values were nonetheless low (0.057 for TRPC1 and 0.076 for S100A1) and unlikely to reflect positive selection.

Discussion

This study is the first to investigate gene regulation in a hibernating primate, and also the first to longitudinally sample WAT from the same individuals with collection points chosen to correlate with the most drastic changes in physiology throughout the circannual cycle. Thus, our study succeeds in providing novel insights into this remarkable phenotypic phenomenon. The dynamic nature of the dwarf lemur WAT transcriptome provides an expansive view of the molecular mechanisms associated with the physiological changes that allow a hibernating primate to withstand long periods of physical inertia and a total lack of food intake by relying solely upon stored fat as fuel. Genes involved in metabolic processes, basic cellular processes (e.g., DNA replication and repair, protein trafficking, cell proliferation), cellular adhesion, blood coagulation, and immune response were found to show highly variable expression profiles that correlate with physiological state (Torpor vs. Active). However, we acknowledge that based on lack of significant enrichment from our functional enrichment analysis, our results should be interpreted carefully. Based on our stringent cut-offs in determining differential expression and, therefore, relatively few numbers of DE genes, we cannot rule out that other genes in the same functional groupings do not show opposite patterns.

Hibernating mammals switch from carbohydrates to lipids as their primary fuel source during the winter hibernation season when seasonal food resources become scarce (Boyer and Barnes 1999; Buck et al. 2002; Carey et al. 2003). Consistent with this physiological transition, we find significant seasonal changes of metabolic gene expression in WAT, including overexpression of fatty acid catabolic genes and lipid transport genes during the winter fasting period. We find instances of other genes involved in fatty acid oxidation that are significantly upregulated during Torpor relative to the Active state. For example, ACSL1 belongs to a family of enzymes that catalyze the esterification of fatty acids with CoA. It is suggested that members of the ACSL gene family drive the uptake of long-chain fatty acids for β-oxidation, and are required for cold thermogenesis (Ellis et al. 2010). Additionally, AASS, a gene product that catalyzes the first two steps in the lysine degradation pathway (Tondo et al. 2013; Sacksteder et al. 2000), is upregulated during Torpor in dwarf lemur WAT. This is important as a significant amount of metabolic energy can come from amino acid metabolism, especially during starvation conditions, as experienced by hibernating animals. Conversely, we uncover expression profiles of genes involved in adipogenesis and fatty acid biosynthesis to be significantly downregulated during Torpor. ADIRF, a gene that is thought to be solely expressed in adipose tissue, putatively aids in the activation of transcription factors involved in adipogenesis, and contributes to increased glucose uptake in adipocytes (Ni et al. 2012; Qiu 2013). This gene is overexpressed in human obesity patients (Ni et al. 2012; Qiu 2013), though more functional studies are needed to determine the exact mechanism for its involvement in adipogenesis.

Our study demonstrates that dwarf lemurs show some DE genes that are also differentially expressed in other hibernating species (table 3). It is worth noting, however, that differential expression of specific genes or protein products has historically shown little consistency when comparing data from previous studies (Villanueva-Cañas et al. 2014) and indeed, we find opposite patterns in few of the genes we identified in dwarf lemurs as differentially expressed when comparing them with American black bears and bat species. However, a closer look taking physiology into account shows that some of these results might not be surprising. AASS, involved in lysine degradation, shows differences in regulation patterns when comparing dwarf lemur WAT with the liver of two species of ground squirrel. Liver and WAT have different responses to amino acid degradation, so its unsurprising that they would show variable gene expression patterns during a torpor bout in different tissues surveyed. We also find variable regulatory patterns of PFKP, a gene involved in glycolysis. Lei et al. (2014) find this gene to be downregulated in the brain of torpid bats, in contrast to what we find in dwarf lemurs. This can be partly explained by the fact that brain and adipose tissue have unique metabolic profiles, especially during periods of starvation, as in torpid animals. Under starvation conditions, the brain switches from glucose to ketone bodies as fuel, a notion that is supported by Lei et al’s finding. This is not the case in WAT, as glucose conservation may be less important.

These results speak to the challenging nature of studies that examine patterns of comparative gene regulation, as well as the need for consistent sampling strategies. However, it is well established that genes function through pathways. Previous work from our group suggests that by expanding differential gene expression profiles to include genes belonging to the same molecular pathway, a higher proportion of genes are identified as shared between hibernators (e.g., seeing high concordance among data sets in beta-oxidation pathways) (Villanueva-Cañas et al. 2014). This lends support to the hypothesis that hibernating species across mammalian phylogeny have evolved similar approaches of gene regulation to combat periods of resource scarcity, although further investigations are required to fully elucidate this question. In this line, we also tested the hypothesis that hibernation-related genes could be evolving under a model of positive selection in hibernating species in comparison with other mammals. For that, we estimated the nonsynonymous to synonymous substitution rate ratio (ω) in different branches of a tree containing hibernating and nonhibernating species (supplementary fig. S3, Supplementary Material online). We found that the model with significantly higher ω in hibernating species was only supported for 2 of 12 genes analyzed, which was not significantly different from what would be expected by chance (supplementary file S2, Supplementary Material online). This is in agreement with previous work wherein we specifically tested the hypothesis that hibernation-related genes were enriched in signals of positive selection (Villanueva-Cañas et al. 2014), finding no support for differential signals of positive selection. Both analyses point in the direction of gene regulation as a key player in the hibernation phenotype rather than the modification of existing genes.

Establishing Methodological Approaches for Investigations with Free-Ranging Populations

Captive conditions offer a unique opportunity for systematically collecting tissue samples under carefully controlled settings. Though the complete longitudinally sampled group is small (N = 4 individuals), the DLC holds the world’s only colony of captive dwarf lemurs, thus providing an exclusive resource for these studies and thus yielding the means to use a sampling approach that has, to our knowledge, not yet been used for transcriptomic-based studies on the hibernation phenotype. The observed duration of torpor bouts and reduction in metabolic rates and Tb are not as extreme in captive animals as those found existing under natural conditions (Dausmann et al. 2005, 2009). Numerous studies in other species have suggested that hibernation patterns differ between intraspecific captive and free-ranging populations (Geiser and Ferguson 2001; Geiser et al. 2007). Although hibernation patterns from captive studies are often used to predict models of energetics and survival in wild populations (Bonaccorso and McNab 1997; Bartels et al. 1998; Armitage et al. 2003; Geiser et al. 2007; Karpovich et al. 2009), the physiology and behavior of hibernation are strongly affected by acclimation to captivity. Nonetheless, captive individuals can be induced to enter bouts of torpor during the winter months through photoperiodic and temperature alterations, with sampling strategies carefully planned for capturing the critical physiological states of interest. Therefore, data from captive studies should be extrapolated to free-ranging populations with caution.

Replicate studies of free-ranging populations in their natural habitat will be crucial for fully understanding the underling genetic mechanisms of primate hibernation. Additionally, leveraging comparative studies of gene expression data from several mammalian hibernators can provide important clues into the specific genetic controls associated with the hibernation phenotype, in addition to providing insight to questions associated with the evolution of the hibernation phenotype within Class Mammalia.

Ethics Statement

Animal handling was carried out in strict accordance with the approval of Duke University’s Institutional Animal Care and Use Committee (IACUC protocol #A274-10-10).

Data Accessibility

Raw sequence data were deposited into the NCBI Short Read Archive with accession number SRP049327. Additional data were deposited into Dryad (doi:10.5061/dryad.057b7; http://datadryad.org/review?doi=doi:10.5061/dryad.057b7, last accessed 13 Jul 2016) and nucleotide alignments for positive selection analysis are uploaded to Figshare (doi: 10.6084/m9.figshare.2013138; https://figshare.com/s/01b4838218146c6187e5, last accessed 13 Jul 2016).

Supplementary Material

Supplementary files S1 and S2, tables S1–S3, and figures S1–S3 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).