Skip to main content

Gut transcriptomic changes during hibernation in the greater horseshoe bat (Rhinolophus ferrumequinum)

Abstract

Background

The gut is the major organ for nutrient absorption and immune response in the body of animals. Although effects of fasting on the gut functions have been extensively studied in model animals (e.g. mice), little is known about the response of the gut to fasting in a natural condition (e.g. hibernation). During hibernation, animals endure the long term of fasting and hypothermia.

Results

Here we generated the first gut transcriptome in a wild hibernating bat (Rhinolophus ferrumequinum). We identified 1614 differentially expressed genes (DEGs) during four physiological states (Torpor, Arousal, Winter Active and Summer Active). Gene co-expression network analysis assigns 926 DEGs into six modules associated with Torpor and Arousal. Our results reveal that in response to the stress of luminal nutrient deficiency during hibernation, the gut helps to reduce food intake by overexpressing genes (e.g. CCK and GPR17) that regulate the sensitivity to insulin and leptin. At the same time, the gut contributes energy supply by overexpressing genes that increase capacity for ketogenesis (HMGCS2) and selective autophagy (TEX264). Furthermore, we identified separate sets of multiple DEGs upregulated in Torpor and Arousal whose functions are involved in innate immunity.

Conclusion

This is the first gut transcriptome of a hibernating mammal. Our study identified candidate genes associated with regulation of food intake and enhance of innate immunity in the gut during hibernation. By comparing with previous studies, we found that two DEGs (CPE and HSPA8) were also significantly elevated during torpor in liver and brain of R. ferrumequinum and several DEGs (e.g. TXNIP and PDK1/4) were commonly upregulated during torpor in multiple tissues of different mammals. Our results support that shared expression changes may underlie the hibernation phenotype by most mammals.

Background

Some mammals have evolved hibernation to cope with limited food availability and high metabolic energy demands in winter. During hibernation, hibernators spend most of their time in a state known as torpor characterized by dramatic reductions in metabolism, body temperature, heart rate, and oxygen consumption [1, 2]. Periods of torpor are interspersed by metabolically active rewarming phase known as interbout arousals [1] which occur in almost all hibernating mammals [3]. Knowledge about functional changes of organs for adaptation to long-term fasting and hypothermia may provide novel insights into human medicine (e.g. organ preservation).

The gut is the major organ for food digestion and nutrient absorption in the body of animals. Comparing with other organs, the gut is affected greatly by food deprivation during hibernation [4]. It has been shown that fasting can cause atrophy of the intestine with reduction in villus height and crypt depth ( e.g. ground squirrels, [5]; bats, [6]). Studies about the effects of food deprivation and internal food stimuli on the gut have identified important genes, proteins, and molecules which modulate food intake (reviewed in [7], see also [8]). However, most of these studies focused on animals under daily food intake fluctuations. Up to now, studies on the effects of a prolonged fast on the gut function have been rarely conducted in a naturally occurring state, such as hibernation.

Another primary gut function is immunity. It is known that the gut is the largest immune organ in the body of animals [9]. It works the first defense to prevent exogenous pathogens from entering the host cells and tissues [10]. A prolonged fast (e.g. 5–8 months for ground squirrels, [11]) and hypothermia during hibernation have dramatic effects on the intestinal immunity [12]. Although previous studies revealed a substantial decline in the innate and adaptive immune function in multiple tissues during hibernation [reviewed in 13], the classical complement pathway was still active in the intestine of thirteen-lined ground squirrels [13]. Also in ground squirrels, a significant increase of intestinal serum albumin (ALB) expression was observed in winter arousal comparing with summer active [14]. ALB has antioxidant and anti-inflammatory properties [15] and may be involved in innate defense [16]. Previous studies have also shown that interbout arousals from torpor can initiate immune response in ground squirrels [17] and can enhance immune responses to a fungal pathogen in bats [18]. However, little is known about the effect of arousals during hibernation on the intestine immune system. Consequently, more comparative studies will be needed to understand the molecular mechanism of how the intestinal immune system responds to hibernation, in particular to the repeated torpor-arousal cycles during hibernation.

Using proteomic analysis, several differentially expressed proteins have been identified to characterize the hibernating phenotype in the intestine of ground squirrels [14]. However, this previous study only focused on differentially expressed proteins between interbout arousal and summer active. So far, little has been known about genes or proteins whose functions are important in the intestine during torpor. In addition, the previous proteomic screen analysis was only limited to a small number of most abundant and soluble proteins [14]. Thus, transcriptome of intestine based on expression of thousands of genes is quite necessary to investigate molecular interactions underlying the hibernation phenotype in this organ. Recently, comparative transcriptomics have been widely used to identify differentially expressed genes between torpid and active mammals from different species and tissues, such as ground squirrels [19,20,21], bats [18, 22,23,24], primates [25, 26], marsupials [27], and bears [28, 29]. In this study, we present the first case of comparative transcriptomics on the gut in a hibernating species of bats.

Almost all insectivorous bats inhabiting in northern latitudes use hibernation to survive the winter. Some of them hibernate for up to eight months without food supply and their body temperatures can be reduced to be below 10 °C [30]. The greater horseshoe bat (Rhinolophus ferrumequinum) has a wide distribution across Europe, Africa and Asia [31] and it has become a model system in studying hibernation of bats. Up to now, comparative transcriptomic studies have been conducted on the brain and liver between winter torpid and summer active bats of R. ferrumequinum [22, 23]. Additionally, a recent study on this species has indicated that physiological changes during hibernation might alter the gut microbiota community [32]. However, little is known about the response of the gut to hibernation in this species. Here we use comparative transcriptomics to investigate functional changes of the gut during hibernation. In particular, we focus on two primary gut functions: food intake regulation and immunity. By examining the expression differences of thousands of genes among bats from winter torpid, arousal and summer active states, we test three hypotheses: 1) In response to the stress of luminal nutrient deficiency during hibernation, genes and/or gene networks related to regulation of food intake would be upregulated in winter torpor relative to summer active; 2) In order to protect the gut from microbial and inflammatory damage, genes involved in innate immune system will be upregulated in winter torpor and arousal of bats; 3) In order to recover from the torpid state quickly, genes associated with cell proliferation (e.g. cell cycle and cell division) would be overexpressed in winter active bats relative to torpid bats. Moreover, we compared our current results from the gut with previous findings in the brain and liver of R. ferrumequinum to investigate molecular effects of torpor on different tissues of hibernating animals. We predicted that differentially expressed genes involved in coping with hypothermia stress would be shared by all of the three tissues.

Methods

Experimental design and sampling

In this study we captured 15 adult male individuals of the greater horseshoe bat (Rhinolophus ferrumequinum) from Fangshan Cave (39°48′ N, 115°42′ E) in Beijing, China. Three of them were sampled in September when they were active (called Summer Active bats). These active bats were captured when they fly out of the cave just before nightfall. The rest of twelve bats were collected in February when they were torpid with rectal temperatures ranging from 8 °C to 10 °C. Rectal temperature was measured by thermometer (Center 309 data logger). Six of these torpid bats were immediately sacrificed on site (called Torpid bats) when their rectal temperatures were still below 10 °C. To investigate the effects of interbout arousals on animals during hibernation, we imitated the arousal state in the wild by gradually allowing torpid bats to restore their euthermic temperature. Specifically, it took us about two hours to transport the six torpid bats to a hotel near the sampling site. Then we kept these bats in a 27 °C room with air condition and did not supply them any food and water. Three of them were sacrificed in about four hours after capturing when their rectal temperatures were around 20 °C (called Arousal bats). The remaining torpid bats were sacrificed in 24 h after capturing when their rectal temperatures were around 37 °C (called Winter Active bats). Thus, in this study we sampled bats from four different physiological states including three winter states (Torpor, Arousal and Winter Active) and one summer active state (Fig. 1a). It is notable that Arousal and Winter Active here are not truly physiological states occurring in the wild.

Fig. 1
figure1

Experimental design and clustering of the 15 R. ferrumequinum individuals based on expression data of 14,009 genes. a Four physiological states included in this study (Torpor, Arousal, Winter Active, and Summer Active) with their corresponding rectal temperatures; b Principal component analysis (PCA) showing gene expression distance of all individuals; c Hierarchical clustering analysis based on pairwise Pearson correlation of gene expression in all individuals

For tissue sampling, bats were rapidly sacrificed by cervical dislocation. Then the small intestine near the stomach was taken and transferred to RNase-free PCR tubes. To reduce the amount of gut microbiota included in the host tissue, we collected a small piece of intestine tissue (~ 1 cm). In addition, to reduce the effect of tissue heterogeneity on estimates of gene expression, for each individual we collected two samples from different section of the small intestine. Thus, a total of 30 tissue samples were included in this study (15 individuals × 2 replicates). Tissues were fresh frozen in liquid nitrogen and stored in a − 80 °C freezer.

RNA extraction, library construction, sequencing and reads trimming

Total RNA was extracted from each of the 30 samples using TRIzol (Life Technologies Corp., Carlsbad, CA, USA) under the manufacturer’s protocol. cDNA library constructions and sequencing on the Illumina Hiseq X ten (paired-end 150 bp) were conducted in Shanghai Majorbio Bio-Pharm Technology Co. Ltd. (four samples) and OE BioTech (26 samples). Raw sequencing reads were trimmed with TRIMMOMATIC version 0.36 [33] using a sliding window of 4 bp with minimum average PHRED quality score of 20 and minimum reads length of 50 bp. A total of 463 Gb clean data were obtained from all 30 samples with an average of 48.9 million clean reads per sample (Table S1). All filtered reads have been deposited into the NCBI Sequence Read Archive database (SRA) (BioProject ID: PRJNA638701).

Reference transcriptome assembly

We used multiple procedures to obtain a high quality reference transcriptome. First, we generated five assemblies based on 12 libraries from Torpid bats, six from bats in each of the other three states (Arousal, Winter Active and Summer Active), and the combination of all 30 libraries, respectively, using Trinity version 2.8.4 [34] with default settings. For each assembly, redundancy was removed using CD-HIT-EST [35] with 95% sequence similarity threshold and the longest open reading frames (ORFs) were extracted using TransDecoder-5.0.0 (http://transdecoder.github.io). Second, transcripts from each of the five assemblies were combined and redundancy was removed again. Third, the resulting assembly was functionally annotated by searching against the UniProtKB/Swiss-Prot and Pfam database (accessed 16th February 2019) using BLASTx and BLASTp with an E-value cutoff of 1.0E-6. Annotation was further performed by searching against the NCBI non-redundant protein database (nr) database (accessed 16th February 2019) using Diamond version 0.9.24 [36] with a set E-value of 1.0E-6. To eliminate exogenous transcripts, the best hits from the two databases above were used to assign taxonomy to each transcript based on NCBI’s prot.accession2taxid and taxdump mapping file (ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/) using R package taxonomizr. Gene Ontology (GO) terms for each transcript were retrieved and used for downstream functional enrichment analysis. Only transcripts both functionally annotated and assigned to mammals were included in the final reference transcriptome which comprises of 14,009 genes (22,847 transcripts) with 96.4% of them (13,502) assigned with a GO term (Table S2). Completeness of this reference transcriptome was assessed with BUSCO (Benchmarking Universal Single-Copy Orthologs, [37]) by searching against 4104 conserved orthologous genes from 50 mammalian species. A total of 91.10% BUSCOs were detected in the final reference transcriptome (Table S2).

Differential expression analysis

Clean reads from each sample were aligned to the reference transcriptome using BOWTIE2 [38] and abundance of transcripts was estimated using RSEM [39] with default settings. Raw read counts across all samples were combined into an expression matrix. Because all samples were sequenced on two different lanes, sequencing batch effect was checked and adjusted using the R package SVA [40]. Six surrogate variables were evaluated by comparing a full model matrix (~ four stages) with a null model matrix (~ intercept term). Prior to differential expression (DE) analysis, we filtered those low expressed transcripts with the average read depth < 10. Then we used principal component analysis (PCA) and hierarchical clustering analysis to visualize gene expression variance across all 30 samples based on SVA-adjusted matrix. A pre-PCA analysis revealed that samples from each of the four states (Torpor, Arousal, Winter Active, and Summer Active) grouped together (Figure S1). Despite the existence of the tissue heterogeneity across samples (Figure S1), the difference of the two sampling replicates from the same individual was not significant (Student’s t test, all P > 0.05), which indicated that variances of gene expression across all samples were mainly caused by different physiological states. Then we combined the aligned bam files of the two samples from the same individual and quantified the transcripts abundance of the specific individual. Using this method, we re-evaluated expression variance across all 15 individuals using PCA and hierarchical clustering analysis as above.

The inclusion of the two non-independent replicates from the same individual during DE analysis may cause the pseudoreplication problem and underestimate the within-group variation of individuals in each physiological state. In order to reduce the type I error rate caused by pseudoreplication problem [41], we account for the correlation of variation within individual using a randomized blocking design. Specifically, we randomly assign the two samples from the same individual to “PartA” and “PartB” and incorporate the blocking design into differential expression tools with ~sampled part + physiological states. By fitting design model above, we can compare the four physiological states with adjustment of the baseline differences from tissue heterogeneity.

Two methods, edgeR [42] and DESeq2 [43], were used to identify differentially expressed genes (DEGs). To reduce false positives, a gene was considered to be differentially expressed only when it was significant in both edgeR and DESeq2 with a false discovery rate (FDR) < 0.05 and log2 transformed fold change (Log2FC) of > 1 or < − 1. To determine the biological role of DEGs identified in each comparison, we performed functional enrichment analysis using the REACTOME, KEGG pathway database and Gene Ontology (GO) [44] in the R/BIOCONDUCTOR package clusterProfiler [45]. Significance is determined with the Bonferroni corrected P value of < 0.01.

Weighted gene co-expression network analysis (WGCNA)

We next identified modules (groups of co-expressed genes across samples) that differentiate each of the four sates in an unbiased way using a weighted gene co-expression network analysis (WGCNA [46];). Our WGCNA focused on 1614 DEGs identified in all pairwise comparisons of the four states above (see Results). Specifically, the correlation between genes’ expression was transformed to connection strengths using a power function with a thresholding power of β = 12. The resulted adjacency matrixes were used to calculate topological overlap matrix (TOM) and corresponding dissimilarity matrix (1-TOM). Then a cluster dendrogram was built by a dynamic tree cutting algorithm [47] to identify modules with a minimum 25 genes per module. Similar modules were merged using a module dissimilarity threshold of 0.25. The module eigengene [48] was used to summarize the expression profiler of a given module. Association between modules and traits were quantified by correlate eigengenes (summary profiles for each module) and evaluated by BH adjusted P-values (q-value, < 0.01). We calculated the module membership value of a gene (kME, the Pearson correlation between its expression value and the module eigengene) which was used to identify the most important intra-module hub genes associated with phenotypic traits [47, 48]. To determine the biological functions of genes in each module, we performed functional enrichment analysis using the same procedures as above. To further investigate the molecular interactions of genes in each module, we conducted protein-protein interaction (PPI) network analysis with STRING version 11.0 [49] and cytoscape version 3.8.0. (https://cytoscape.org/).

Results

Differential expression analysis

The principal component (PCA) and clustering analysis both revealed that individuals from each of the four states (Torpor, Arousal, Winter Active, and Summer Active) grouped together respectively (Fig. 1b and c). To account for gene expression variation of the two replicates from the same individual (Figure S1), we included all 30 samples in the following DE analysis. However, to reduce the type I error rate caused by pseudoreplication problem, we used a randomized blocking design when performing DE analysis (see details in Materials and methods). Under this blocking design, EdgeR and DESeq2 identified almost the same number of differentially expressed genes (DEGs) in each comparison (Figure S2, Table S3 and S4). In this study only shared DEGs identified by both methods were considered as true DEGs (Fig. 2a and b). After removing the overlapped DEGs identified in all six comparisons, a total of 1614 DEGs were retained (Fig. 2b, Table S5).

Fig. 2
figure2

a Bar graph showing the number of differentially expressed genes (DEGs) identified in each pairwise comparison of the four states. b Venn diagram showing the number of DEGs shared between each pairwise comparison of four physiological states. TP vs SAT: Torpor vs Summer Active; WAT vs SAT: Winter Active vs Summer Active; AS vs SAT: Arousal vs Summer Active; WAT vs TP: Winter Active vs Torpor; AS vs TP: Arousal vs Torpor; WAT vs AS: Winter Active vs Arousal

For upregulated DEGs, significant enrichments of GO terms and/or KEGG and Reactome pathways were detected in Torpor, Winter Active and Arousal relative to other states (Figure S3, Table S6 and S7). We found that those DEGs upregulated in Torpor relative to Summer Active were associated with insulin-like Growth Factor (IGF), and those DEGs upregulated in Winter Active relative to Arousal were associated with carbohydrate metabolism (Figure S3 and Table S7). Those DEGs upregulated in Arousal and Winter Active comparing with Torpor or Summer Active were associated with metabolism of lipids, carbohydrate, cofactors, and vitamins, biological oxidations, digestion, apical plasma membrane, and carbohydrate transmembrane transporter activity (Figure S3, Table S6 and S7). For downregulated DEGs, we found significant enrichments of GO terms (Table S6) and/or KEGG and Reactome pathway (Figure S4 and Table S7) in five comparisons. Among those downregulated in Torpor, Arousal and Winter Active relative to Summer Active (n = 369), a large proportion of them (230, 62.3%) were associated with immune function (Table S6 and S7). Those downregulated in Arousal and Winter Active relative to Torpor were associated with metabolism of proteins, hemostasis and complement and coagulation cascades (Figure S4, Table S6 and S7).

WGCNA

WGCNA based on 1614 DEGs yielded ten modules with the number of genes ranging from 72 to 301 per module (Fig. 3a, b; Table S5). Among these modules, six were found to be positively associated with three winter states (TP-M9 and TP-M10 modules for Torpor; AS-M1, AS-M2 and AS-M3 modules for Arousal; WAT-M4 module for Winter Active). For Summer Active state, we identified four significantly positive modules (SAT-M5, SAT-M6, SAT-M7, and SAT-M8 modules).

Fig. 3
figure3

Weighted Gene Coexpression Network Analysis (WGCNA). a Cluster dendrogram built based on 1614 DEGs in WGCNA. All modules and merged modules with a module dissimilarity threshold of 0.25 were shown in upper and down panel, respectively. b Heatmap showing correlations between module eigengenes and each of the four states. The color scale indicates the strength of correlation. Number of genes included in each module was shown in the bar plot right; c-e Results of protein-protein interaction (PPI) network analysis on DEGs in three WGCNA modules (TP-M9, TP-M10 and AS-M1 modules; c, d and e). Size of nodes indicates module membership values (MM) of specific DEGs in corresponding WGCNA modules. Size of the label indicates the mean expression level (log2CPM) across the samples. Proteins discussed in Results are cycled with red color. Edge transparency indicates the interaction score between two proteins (less transparent stronger correlation of the two proteins)

We found significant GO terms and/or pathways in genes from three modules (Figure S5, Table S8). Specifically, genes in TP-M9 and/or TP-M10 module were associated with Insulin−like Growth Factor (IGF), protein activation cascade, and hemostasis, whereas genes in AS-M1 module were enriched in the process of cell cycle and mitosis. We further conducted PPI network analysis on genes from modules associated with the three winter states to investigate molecular interactions of these genes. PPI networks were identified among genes from three modules, which will be discussed in details below.

First, one network including 19 proteins was constructed in TP-M9 module (Fig. 3c and Table S9). Among these proteins, eight of the top 10 proteins with the largest expression difference between Torpor and Summer Active belong to acute-phase proteins (APPs) (AHSG, ALB, FETUB, FGA, FGB, HP, HPX, and TF). Another two are heat shock proteins (HSPA1A and HSPB1). All of these APPs are in the list of the top 10 proteins with the largest WGCNA module membership (MM) value. The rest two, CPS1 and ITIH3, are associated with urea cycle and hemostasis, respectively. In addition, seven APPs are also in the list of the top 10 proteins with the largest degree of connection (K) value in PPI network. The rest three, CST3, CSF1 and A2M, are associate with protective function and hemostasis. Last, seven proteins above (A2M, ALB, CST3, CSF1, HSPA1A, HSPB1, and TF) are also in the list of the top 10 proteins with the highest mean expression value across Torpid samples.

Second, two networks including 66 proteins were generated in TP-M10 module (Fig. 3d and Table S10). For Network-1 (48 proteins), we mainly focused on those proteins in the list of the top 15 with either the largest expression difference between Torpor and Summer Active, the highest mean expression value across Torpid samples, the largest MM value or K value. Thirty-six such proteins are included. Among them, five and two proteins are associated with regulation of food intake (CCK, CPE, DRD2, GPR17, and VGF) and digestion (SSTR2 and APOA1), respectively. All of these seven proteins are centered by KNG1 which shows the largest K value (Fig. 3d, Network-1). TEX264 is also centered by KNG1 and related to ER-phagy. One (CYR61) may be involved in restoring mucosal homeostasis after intestinal injury [50]. Three (MUC1, AGR2 and TFF3) are related to maintenance and repair of the intestinal mucosa. AGR2 and TFF3 also have the largest mean expression value across Torpid samples. It is notable that another intestinal defense protein (MUC2) and a histone protein (HIST2H2AA3) show the largest expression difference between Torpor and Summer Active, and the later histone protein also exhibits the largest MM value. Another one protein (EDN1) has an effect on ischemic stroke [51]. We also found two insulin-like growth factors (IGF2 and IGFBP7). For Network-2, a heat shock protein (HSPA8) and an iron storage protein ferritin (FTH1), have the largest mean expression value across Torpid samples, and TXNIP, another highly expressed proteins (top 5 in the list of 66 proteins above), but not in the network, is a potent antioxidant.

Third, we identified one network in AS-M1 module including 41 proteins (Fig. 3e and Table S11). Among the top 10 upregulated genes in Arousal vs. Torpor, nine are involved in cell cycle and mitosis (AURKB, BIRC5, BUB1B, CDCA3, CCNB1, NUSAP1, RRM2, TOP2A, and ZBTB16). The left one, KIF19 plays an important role in regulation of the length of motile cilia. NFKBIA has the largest MM value and is related to immune function. Among the other top 10 proteins with the largest MM, five were also in the list of the top 10 upregulated genes in Arousal vs. Torpor (AURKB, BIRC5, CCNB1, NUSAP1, and RRM2), and together with another four proteins (CCNF, KIF4A, MYBL2, and SLC25A5), are all involved in cell cycle and mitotic process. The left one, KIF19 was also in the top 10 upregulated genes in Arousal vs. Torpor as described above. Among the above proteins, SLC25A5 has the highest expression value in Arousal. Last, all of the top 10 proteins with the largest K value are associated with cell cycle and mitosis.

Discussion

To our knowledge, this is the first study to assess effects of hibernation on the gut of animals in natural conditions using comparative transcriptomics. We identified 1614 differentially expressed genes (DEGs) during four physiological states (Torpor, Arousal, Winter Active, and Summer Active). Using a weighted gene co-expression network analysis (WGCNA), we assigned 926 DEGs into six modules associated with the three winter states. Our current results include all of the proteins (ALB, HSPA8, APOA1, and HMGCS2) that showed over two-fold expression changes in winter Arousal vs. summer Active in the intestine of ground squirrel [14]. This consistency supports the reliability of our current results.

In this study we observe overexpression of genes in winter Torpor or Arousal whose functions are related to several important aspects of the gut, such as intestinal defense (MUC1, MUC2, AGR2, and TFF3), cilia formation (KIF19), responses to ischemia (EDN1), and cell cycle and mitosis (see details in Results). In order not to be lost in the map, here we focus on transcriptomic changes of two primary gut functions: regulation of food intake and immune response.

Regulation of food intake in the gut during torpor

Because the gut is the major organ for nutrient absorption in animals, we predict that the deficiency of luminal nutrient during hibernation may induce overexpression of those genes that reduce food intake or suppress the appetite. Consistent with this prediction, we find that at least six genes in TP-M10 module (CCK, CPE, GPR17, VGF, RXFP4, and DRD2) are involved in regulation of food intake. Proteins encoded by these genes are mostly centered by KNG1 in the protein-protein interaction network (PPI network, Fig. 3d and they were all upregulated at least four-fold in Torpor vs. Active. Below we will discuss their functional roles in food intake regulation one by one.

CCK (cholecystokinin) is one of the well-known gut hormones which act as important satiety signals mediator in controlling food intake and regulating absorption of nutrients [52, 53]. It is known that the release of CCK is stimulated by nutrient signals from food digestion [54, 55]. However, these nutrient signals may be blocked during the long-term fasting in hibernating animals. Recently, gastrointestinal (GI) distension, as an afferent signal, has been proved to be more sensitive in regulating food intake than nutrient-stimulated signals in GI tract [56]. Thus, GI distension may serve as a mediator transferring satiety signals and stimulate the release of CCK [56, 57]. Interestingly, CCK is the only upregulated gene among the three hormones secreted from the intestine [CCK, GIP and PYY, 59] that induce satiety signals to the nucleus tractus solitarius (NTS) though vagal afferents [58]. Among them, GIP was downregulated in Torpor relative to Active, and PYY was not differentially expressed in this study. In addition, CCK can modulate the sensitivity of gastric afferents to leptin [59] and inhibit gastric emptying in animals [60]. In rats CCK mediates satiety signals by inducing the release of leptin stored in stomach [61]. Leptin, a steroid hormone, in return, can enhance satiating effect of CCK [57] and the role of leptin in decreasing food intake has been well documented in mammals [8]. Later studies have proposed that there is a synergistic interaction between CCK and leptin to reduce food intake [7, 62]. Consistent with this proposal, we found that LEP, encoding leptin, was expressed with the second largest fold change among all upregulated genes in Torpor relative to Active (nearly 20-fold) although it was identified as the DEG only by DESeq2 (Table S4). Consistent with our finding, CCK was also significantly elevated in Torpor vs. Active in the liver of R. ferrumequinum ([23], Table S12]). CPE (carboxypeptidase E), an obesity susceptibility gene, has been shown to be responsible for the processing of pro-CCK into its bioactive forms [63]. Stimulation of CPE can increase energy efficiency [64] which is essential when energy is limited during hibernation. The importance of CPE for hibernating animals is further supported by the observation that it was upregulated in Torpor vs. Active in brain and liver of R. ferrumequinum ([22, 23]; Table S12). GPR17, a G protein-coupled receptor, acts as a target of FoxO1 to regulate appetite and energy expenditure [65]. The expression of FOXO1 can be induced in response to fasting [66, 67]. However, in the current study FOXO1 was not differentially expressed in all six comparisons and it maintained a high expression level in each of the four physiological states with a mean log2CPM value of 5. STAT3, an antagonism to FoxO1, also maintained a high expression level in all of the four states and showed the highest expression in Torpor (mean log2CPM > 7). Leptin can inactivate FoxO1 and reduce expression of GPR17 in nucleus. The inactivation of FoxO1 allows STAT3 to inhibit the expression of AGRP and stimulate the expression of POMC [68], which results in increased sensitivity to insulin and leptin and reduced food intake [68]. Additionally, FoxO1 also regulates CPE to integrate food intake with energy expenditure [64]. In summary, our current results suggest that long-term fasting and energy depression during hibernation may first induce the release of endogenous hormones (e.g. CCK and Leptin) which can help to regulate food intake by triggering the circuit between central nerve system and the intestinal afferent vagus nerve. To clearly illustrate the role of genes/proteins discussed above in food intake regulation, we summarized their functional interactions as in Fig. 4a (Modified from [57, 68]).

Fig. 4
figure4

Schematic view of the molecular interactions of genes/proteins involved in response to the stress of luminal nutrient deficiency during hibernation. a The molecular interactions of DEGs in energy supply and regulation of food intake by transferring satiety signals between central nerve system and the intestinal afferent vagus nerve, and Leptin signaling in the ARC (arcuate nucleus). b Cluster of six DEGs (CCK, CPE, DRD2, GPR17, TEX264, and VGF) associated with food intake regulation based on expression patterns across the four states

VGF, a neuronal and endocrine peptide precursor, and its derived peptides can decrease food intake [69,70,71]. It was notable that this gene was significantly downregulated in Torpor vs. Active in the brain of R. ferrumequinum [22]. This may be caused by differential roles of VGF in the gut and brain [72]. Another G protein-coupled receptor, RXFP4, is expressed in the enteric nervous system and its ligand is the gut novel hormone insulin-like peptide 5 (INSL5). A recent study found that fasting and long-term energy deprivation can induce the expression of INSL5 which enhances appetite by stimulating RXFP4 [73]. It was notable that INSL5 was almost not expressed in the gut. The last one, DRD2 can also mediate food intake via dopamine signaling [74].

Another way to respond to stress of nutrient deficiency or starvation is to activate endoplasmic reticulum-phagy (ER-phagy). This selective autophagy helps mammalian cell to acquire amino acids and other cellular components by selectively transporting cytosolic contents to lysosomes for degradation. Two recent studies have proved that TEX264 is a critical determinant of ER-phagy in response to nutrient deprivation ([75, 76]; see also Fig. 4). In addition, fasting induces an overexpression of HMGCS2 in Torpor (almost 5-fold in Torpor vs. Active, Table S3 and S4). The encoded protein by this gene is a key enzyme regulating ketone body production. The increased ability for ketogenesis can produce adequate energy supply for the gut. Consistent with our result here, HMGCS2 was also elevated 2.2-fold in winter Arousal vs. summer Active in the intestine of ground squirrels [14].

Finally, we found that six DEGs (CCK, CPE, DRD2, GPR17, TEX264, and VGF) showed a decreased expression level across the four states from Torpor to Summer Active (Fig. 4b), corresponding to an increase of energy demand from food intake along the rise of body temperature. This implies that genes regulating food intake and energy supply may be important in triggering cyclic variation of energy demand during hibernation.

Regulation of immune system in the gut during hibernation

Our current study revealed a large proportion of downregulated DEGs in three winter states relative to Active were associated with immunity (Table S6 and S7), supporting a systemic suppression of immune function during hibernation [77]. However, consistent with our second hypothesis, we found that several upregulated genes in Torpor or Arousal were related to innate immunity.

First, among the 30 proteins in the PPI network of TP-M9 module (Fig. 3c and Table S9), eight of the top 10 upregulated proteins in Torpor vs. Summer Active belong to acute-phase proteins which play important roles in the early response to infection or tissue injury [78]. Among them, AHSG (also known as FETUA) and FETUB belong to the set of negative acute-phase genes [79] and they play anti-inflammatory roles by counter-regulating the innate immune response [80]. It is notable that AHSG shows the largest expression difference between Torpor and Summer Active and also the largest K value. ALB has the second largest K value and its protein (Albumin) has antioxidant and anti-inflammatory properties [15]. This protein exhibited the greatest fold change in the intestine of ground squirrels in winter Arousal vs. summer Active [14]. Similarly, APOA1 was also overexpressed in hibernating ground squirrel intestine [14] and has potential effects of anti-inflammatory. Hemopexin encoded by HPX exhibits the biggest MM value. This protein is a heme-binding plasma glycoprotein and can act as an antioxidant against oxidative damage mediated by hemoglobin [81]. Haptoglobin encoded by HP also acts as an antioxidant and anti-inflammatory [82]. Similar to hemopexin, it also plays a crucial role in against hemoglobin-mediated oxidative damage [83]. Another important protein in the network of TP-M9 module is CSF1, macrophage colony-stimulating factor, which can restore innate immunity after injury of liver in mice [84]. Last, KNG1 has the largest K value in the networks of proteins in TP-M10 module. This protein has been proved to be involved in inflammation, coagulation, and innate immunity [85]. Consistent with our findings above, previous studies on the intestine and lung of thirteen-lined ground squirrels showed that certain aspects of innate immunity were still maintained during torpor [13, 86]. Another study on the bone marrow of thirteen-lined ground squirrels revealed a shift to innate immune responses during hibernation [21]. Overall, our current and previous studies support that innate immune system plays an essential role in protection against infection or tissue injury during torpor.

Second, among the proteins in AS-M1 module (Table S5), NFKBIA exhibits the largest MM value and this protein has been reported to alleviate gut infection by reinforcing innate immune response via the NF-kB system [87]. Another protein, CEBPD, is an important transcription factor regulating the expression of genes involved in immune and inflammatory responses. This protein has been shown to play an essential role in Gamma-tocotrienol mediated protection against ionizing radiation (IR)-induced hematopoietic and intestinal injury [88]. A further protein (GGT1) also plays an important role in cell antioxidant defense mechanism [89]. Accordingly, our current results support an important role of interbout arousals from torpor in enhancing immune responses during hibernation [17, 18].

It is notable that ACE2, encoded a receptor for human coronaviruses SARS-CoV and 2019-nCoV [90], was downregulated over 2 fold in Torpor relative to Active (Table S4) although currently we have no idea about its effect on the interaction of ACE2 with coronaviruses during hibernation.

Comparing with previous studies

In order to test whether hibernation has a differential effect on different tissues of a hibernating animal, we compared our current results with previous comparative transcriptomics on the brain and liver in R. ferrumequinum [22, 23]. Among upregulated DEGs in Torpor vs. Summer Active, we found two DEGs (CPE and HSPA8) overlapped by the three tissues. In addition, we identified 11, 13 and 32 DEGs shared by the gut and liver, the gut and brain, and the brain and liver, respectively (Table S12).

Consistent with our prediction that shared DEGs by multiple tissues will be important in coping with hypothermia during hibernation, we found one head shock protein (HSPA8) upregulated in all of the three tissues. This heat shock protein has the second largest expression value in TP-M10 module (Table S10, Fig. 3d) and was also significantly overexpressed in winter Arousal vs. summer Active in the intestine of ground squirrels [14]. These suggest that this protein may belong to ‘protective’ protein during hibernation. In addition, other heat shock proteins were found in the list of shared upregulated genes in torpor vs. active between the gut and liver (HSPB1) or between the brain and liver (HSP90B1, HSP90AA1, HSPA4L, and HSPH1) (Table S12), suggesting their general roles in protection against stress during torpor [12].

By comparing with other previous studies, we found that many of the shared upregulated genes in torpor vs. active by two of the three tissues in R. ferrumequinum were also reported to be important during torpor in other hibernating mammals. For example, TXNIP, one of the 11 shared upregulated genes between the gut and liver (Table S12), encodes thioredoxin-interacting protein which is a potent antioxidant. This gene has been shown to be significantly elevated in torpor vs. active in all three tissues (liver, brain and muscle) of small South American marsupials [27] and was also over three fold higher in torpor vs. active in the bone marrow of thirteen-lined ground squirrels [21]. Another example is PDK1, one of the 13 shared upregulated gene between the gut and brain (Table S12). Its encoded protein, PDK1, is a protein kinase. PDK1 and its other isozymes (PDK 2, 3 and 4) phosphorylate pyruvate dehydrogenase (PDH) and inhibits its activity [91, 92], which promotes the shift from carbohydrate metabolism to lipolysis. Suppression of PDH activity by upregulation of PDK isozymes (PDK 1 and 4) has been well documented in multiple tissues of different hibernating mammals (e.g. white adipose tissue in dwarf lemurs, [26], brain in marsupial, [27], liver and skeletal muscle in ground squirrels, [93]). The commonality of overexpression of these genes in torpor in different tissues and different mammals further supports that genes and biochemical process underlying the hibernation phenotype may be shared by all mammals (see also [94]).

Conclusions

This study has three major findings. First, we identified candidate genes that are important in response to the stress of luminal nutrient deficiency during hibernation. Our PPI network results suggest that they are interacted with each other in regulation of food intake. Second, we detected overexpression of genes in Torpor and/or Arousal whose functions are involved in innate immune response. Such genes may be important in preservation of organ function during extended hypothermia. Third, by comparing with previous studies, we find several important genes whose expressions are commonly elevated during torpor in multiple tissues of different hibernating mammals. Thus, our findings in this study can be applicable to other hibernating mammals.

It is now known that the gut microbiota plays essential roles not only in food digestion, energy harvest and metabolism, but also in immune function of their hosts [95, 96]. The host-microbial cross-talk is particularly dynamic in hibernating mammals with a prolonged fast [97]. Previous studies on the gut microbiota in hibernating mammals have shown that the structure and composition of the gut microbiota have changed greatly during hibernation (e.g. in ground squirrels, [98], in horseshoe bats, [32]). In the future, functional comparisons of the structure and composition of the gut microbiota between summer active and winter torpid animals will be important if we aim to understand the role of molecular interactions of the gut and microbiota in response to stressful conditions (e.g. extended fasting and hypothermia) during hibernation.

Availability of data and materials

The raw sequencing reads were deposited in the NCBI Sequence Read Archive (SRA) under Bioproject PRJNA638701, SRR12006225-SRR12006254. Assembly of the reference transcriptome is available in Dryad: https://doi.org/10.5061/dryad.3n5tb2rd8. The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. All data generated or analysed during this study are included in this published article and its supplementary information files.

Abbreviations

BUSCO:

Benchmarking Universal Single-Copy Orthologs

DEGs:

Differentially expressed genes

ER-phagy:

Endoplasmic reticulum-phagy

FDR:

False discovery rate

GO:

Gene Ontology

Log2FC:

Log2 transformed fold change

ORFs:

Open reading frames

MM:

Module membership value

nr:

Non-redundant protein database

PCA:

Principal component analysis

PPI:

Protein-protein interaction

RNA-Seq:

RNA sequencing

SRA:

Sequence read archive database

Tb:

Rectal temperature

TOM:

Topological overlap matrix

WGCNA:

Weighted gene co-expression network analysis

Genes or proteins, ALB:

Albumin

CCK:

Cholecystokinin

HPX:

Hemopexin

HP:

Haptoglobin

CSF1:

Macrophage colony-stimulating factor

PDH:

Phosphorylate pyruvate dehydrogenase

TXNIP:

Thioredoxin-interacting protein

GPR17:

G protein-coupled receptor 17

References

  1. 1.

    Carey HV, Andrews MT, Martin SL. Mammalian hibernation: cellular and molecular responses to depressed metabolism and low temperature. Physiol Rev. 2003;83:1153–81.

    CAS  PubMed  Google Scholar 

  2. 2.

    Geiser F. Metabolic rate and body temperature reduction during hibernation and daily torpor. Annu Rev Physiol. 2004;66:239–74.

    CAS  PubMed  Google Scholar 

  3. 3.

    Van Breukelen F, Martin SL. The hibernation continuum: physiological and molecular aspects of metabolic plasticity in mammals. Physiology. 2015;30:273–81.

    PubMed  Google Scholar 

  4. 4.

    McCue MD. Starvation physiology: reviewing the different strategies animals use to survive a common challenge. Comp Biochem Physiol A Mol Integr Physiol. 2010;156:1–18.

    PubMed  Google Scholar 

  5. 5.

    Carey HV, Cooke HJ. Effect of hibernation and jejunal bypass on mucosal structure and function. Am J Physiol Gastrointest Liver Physiol. 1991;261:G37–44.

    CAS  Google Scholar 

  6. 6.

    Paksuz EP. The effect of hibernation on the morphology and histochemistry of the intestine of the greater mouse-eared bat, Myotis myotis. Acta Histochem. 2014;116:1480–9.

    CAS  PubMed  Google Scholar 

  7. 7.

    Cummings DE, Overduin J. Gastrointestinal regulation of food intake. J Clin Invest. 2007;117:13–23..

    CAS  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Florant GL, Healy JE. The regulation of food intake in mammalian hibernators: a review. J Comp Physiol B. 2012;182:451–67.

    CAS  PubMed  Google Scholar 

  9. 9.

    Kraehenbuhl J-P, Neutra MR. Molecular and cellular basis of immune protection of mucosal surfaces. Physiol Rev. 1992;72:853–79.

    CAS  PubMed  Google Scholar 

  10. 10.

    Mathew AG. Nutritional influences on gut microbiology and enteric diseases. Biotechnology in the Feed Industry. In: Lyons TP, ed. 2001.

  11. 11.

    Kurtz CC, Carey HV. Seasonal changes in the intestinal immune system of hibernating ground squirrels. Dev Comp Immunol. 2007;31:415–28.

    CAS  PubMed  Google Scholar 

  12. 12.

    Carey H, Mangino M, Southard J. Changes in gut function during hibernation: implications for bowel transplantation and surgery. Gut. 2001;49:459–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Maniero GD. Classical pathway serum complement activity throughout various stages of the annual cycle of a mammalian hibernator, the golden-mantled ground squirrel, Spermophilus lateralis. Dev Comp Immunol. 2002;26:563–74.

    CAS  PubMed  Google Scholar 

  14. 14.

    Martin SL, Epperson LE, Rose JC, Kurtz CC, Ané C, Carey HV. Proteomic analysis of the winter-protected phenotype of hibernating ground squirrel intestine. Am J Physiol Regul Integr Comp Physiol. 2008;295:R316–28.

    CAS  PubMed  Google Scholar 

  15. 15.

    Caraceni P, Rosenblum ER, Van Thiel DH, Borle AB. Reoxygenation injury in isolated rat hepatocytes: relation to oxygen free radicals and lipid peroxidation. Am J Physiol Gastrointest Liver Physiol. 1994;266:G799–806.

    CAS  Google Scholar 

  16. 16.

    Shamay A, Homans R, Fuerman Y, Levin I, Barash H, Silanikove N, et al. Expression of albumin in nonhepatic tissues and its synthesis by the bovine mammary gland. J Dairy Sci. 2005;88:569–76.

    CAS  PubMed  Google Scholar 

  17. 17.

    Prendergast BJ, Freeman DA, Zucker I, Nelson RJ. Periodic arousal from hibernation is necessary for initiation of immune responses in ground squirrels. Am J Physiol Regul Integr Comp Physiol. 2002;282:R1054–62.

    CAS  PubMed  Google Scholar 

  18. 18.

    Field KA, Sewall BJ, Prokkola JM, Turner GG, Gagnon MF, Lilley TM, et al. Effect of torpor on host transcriptomic responses to a fungal pathogen in hibernating bats. Mol Ecol. 2018;27:3727–43.

    CAS  Google Scholar 

  19. 19.

    Schwartz C, Hampton M, Andrews MT. Seasonal and regional differences in gene expression in the brain of a hibernating mammal. PLoS One. 2013;8:e58427.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Bogren LK, Grabek KR, Barsh GS, Martin SL. Comparative tissue transcriptomics highlights dynamic differences among tissues but conserved metabolic transcript prioritization in preparation for arousal from torpor. J Comp Physiol B. 2017;187:735–48.

    CAS  PubMed  Google Scholar 

  21. 21.

    Cooper ST, Sell SS, Fahrenkrog M, Wilkinson K, Howard DR, Bergen H, et al. Effects of hibernation on bone marrow transcriptome in thirteen-lined ground squirrels. Physiol Genomics. 2016;48:513–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Lei M, Dong D, Mu S, Pan Y-H, Zhang S. Comparison of brain transcriptome of the greater horseshoe bats (Rhinolophus ferrumequinum) in active and torpid episodes. PLoS One. 2014;9:e107746.

    PubMed  PubMed Central  Google Scholar 

  23. 23.

    Xiao Y, Wu Y, Sun K, Wang H, Zhang B, Song S, et al. Differential expression of hepatic genes of the greater horseshoe bat (Rhinolophus ferrumequinum) between the summer active and winter torpid states. PLoS One. 2015;10:e0145702.

    PubMed  PubMed Central  Google Scholar 

  24. 24.

    Lilley TM, Prokkola JM, Blomberg AS, Paterson S, Johnson JS, Turner GG, et al. Resistance is futile: RNA-sequencing reveals differing responses to bat fungal pathogen in Nearctic Myotis lucifugus and Palearctic Myotis myotis. Oecologia. 2019;191:295–309.

    PubMed  PubMed Central  Google Scholar 

  25. 25.

    Faherty SL, Villanueva-Cañas JL, Blanco MB, Albà MM, Yoder AD. Transcriptomics in the wild: hibernation physiology in free-ranging dwarf lemurs. Mol Ecol. 2018;27:709–22.

    CAS  PubMed  Google Scholar 

  26. 26.

    Faherty SL, Villanueva-Cañas JL, Klopfer PH, Albà MM, Yoder AD. Gene expression profiling in the hibernating primate, Cheirogaleus medius. Genome Biol Evol. 2016;8:2413–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Nespolo RF, Gaitan-Espitia JD, Quintero-Galvis JF, Fernandez FV, Silva AX, Molina C, et al. A functional transcriptomic analysis in the relict marsupial Dromiciops gliroides reveals adaptive regulation of protective functions during hibernation. Mol Ecol. 2018;27:4489–500..

    CAS  PubMed  Google Scholar 

  28. 28.

    Jansen HT, Trojahn S, Saxton MW, Quackenbush CR, Hutzenbiler BDE, Nelson OL, et al. Hibernation induces widespread transcriptional remodeling in metabolic tissues of the grizzly bear. Commun Biol. 2019;2:1–10..

    CAS  Google Scholar 

  29. 29.

    Mugahid DA, Sengul TG, You X, Wang Y, Steil L, Bergmann N, et al. Proteomic and transcriptomic changes in hibernating grizzly bears reveal metabolic and signaling pathways that protect against muscle atrophy. Sci Rep. 2019;9:1–16.

    Google Scholar 

  30. 30.

    Stawski C, Willis C, Geiser F. The importance of temporal heterothermy in bats. J Zool. 2014;292:86–100.

    Google Scholar 

  31. 31.

    Csorba G, Ujhelyi P, Thomas N. Horseshoe bats of the world (Chiroptera: Rhinolophidae). Shrewsbury: Alana books; 2003.

    Google Scholar 

  32. 32.

    Xiao G, Liu S, Xiao Y, Zhu Y, Li A, Li Z, et al. Seasonal changes in gut microbiota diversity and composition in the greater horseshoe bat. Front Microbiol. 2019;10:2247.

    PubMed  PubMed Central  Google Scholar 

  33. 33.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29:644.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.

    CAS  PubMed  Google Scholar 

  36. 36.

    Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59.

    CAS  PubMed  Google Scholar 

  37. 37.

    Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.

    PubMed  Google Scholar 

  38. 38.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28:882–3.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Hurlbert SH. Pseudoreplication and the design of ecological field experiments. Ecol Monogr. 1984;54:187–211.

    Google Scholar 

  42. 42.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    CAS  PubMed  Google Scholar 

  43. 43.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    PubMed  PubMed Central  Google Scholar 

  44. 44.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: J Integrative Biol. 2012;16:284–7.

    CAS  Google Scholar 

  46. 46.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    PubMed  PubMed Central  Google Scholar 

  47. 47.

    Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R. Bioinformatics. 2008;24:719–20.

    CAS  PubMed  Google Scholar 

  48. 48.

    Langfelder P, Mischel PS, Horvath S. When is hub gene selection better than standard meta-analysis? PLoS One. 2013;8:e61505.

  49. 49.

    Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47:D607–13.

    CAS  PubMed  Google Scholar 

  50. 50.

    Choi JS, Kim K-H, Lau LF. The matricellular protein CCN1 promotes mucosal healing in murine colitis through IL-6. Mucosal Immunol. 2015;8:1285–96.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Zhang L, Sui R. Effect of SNP polymorphisms of EDN1, EDNRA, and EDNRB gene on ischemic stroke. Cell Biochem Biophys. 2014;70:233–9.

    CAS  PubMed  Google Scholar 

  52. 52.

    Cordomí A, Fourmy D, Tikhonova IG. Gut hormone GPCRs: structure, function, drug discovery. Curr Opin Pharm. 2016;31:63–7.

    Google Scholar 

  53. 53.

    Crawley JN, Corwin RL. Biological actions of cholecystokinin. Peptides. 1994;15:731–55.

    CAS  PubMed  Google Scholar 

  54. 54.

    Peikin SR. Role of cholecystokinin in the control of food intake. Gastroenterol Clin N Am. 1989;18:757–75.

    CAS  Google Scholar 

  55. 55.

    Moran TH, Schwartz GJ. Neurobiology of cholecystokinin. Crit Rev Neurobiol. 1994;9:1–28.

    CAS  PubMed  Google Scholar 

  56. 56.

    Bai L, Mesgarzadeh S, Ramesh KS, Huey EL, Liu Y, Gray LA, Aitken TJ, Chen Y, Beutler LR, Ahn JS. Genetic identification of vagal sensory neurons that control feeding. Cell. 2019;179:1129–43.

    CAS  PubMed  Google Scholar 

  57. 57.

    Schwartz MW, Woods SC, Porte D, Seeley RJ, Baskin DG. Central nervous system control of food intake. Nature. 2000;404:661–71.

    CAS  PubMed  Google Scholar 

  58. 58.

    Dockray GJ. Enteroendocrine cell signalling via the vagus nerve. Curr Opin Pharm. 2013;13:954–8.

    CAS  Google Scholar 

  59. 59.

    Wang Y, Tache Y, Sheibel AB, Go V, Wei J. Two types of leptin-responsive gastric vagal afferent terminals: an in vitro single-unit study in rats. Am J Physiol Regul Integr Comp Physiol. 1997;273:R833–7.

    CAS  Google Scholar 

  60. 60.

    Scarpignato C, Varga G, Corradi C. Effect of CCK and its antagonists on gastric emptying. J Physiol Paris. 1993;87:291–300.

    CAS  PubMed  Google Scholar 

  61. 61.

    Bado A, Levasseur S, Attoub S, Kermorgant S, Laigneau J-P, Bortoluzzi M-N, Moizo L, Lehy T, Guerre-Millo M, Le Marchand-Brustel Y. The stomach is a source of leptin. Nature. 1998;394:790–3.

    CAS  PubMed  Google Scholar 

  62. 62.

    Barrachina MD, Martinez V, Wang L, Wei JY, Taché Y. Synergistic interaction between leptin and cholecystokinin to reduce short-term food intake in lean mice. Proc Natl Acad Sci. 1997;94:10455–60.

    CAS  PubMed  Google Scholar 

  63. 63.

    Cain BM, Wang W, Beinfeld MC. Cholecystokinin (CCK) levels are greatly reduced in the brains but not the duodenums of Cpefat/Cpefat mice: a regional difference in the involvement of carboxypeptidase E (Cpe) in pro-CCK processing. Endocrinology. 1997;138:4034–7.

    CAS  PubMed  Google Scholar 

  64. 64.

    Plum L, Lin HV, Dutia R, Tanaka J, Aizawa KS, Matsumoto M, et al. The obesity susceptibility gene Cpe links FoxO1 signaling in hypothalamic pro-opiomelanocortin neurons with regulation of food intake. Nat Med. 2009;15:1195.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Ren H, Orozco IJ, Su Y, Suyama S, Gutiérrez-Juárez R, Horvath TL, et al. FoxO1 target Gpr17 activates AgRP neurons to regulate food intake. Cell. 2012;149:1314–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Kim M-S, Pak YK, Jang P-G, Namkoong C, Choi Y-S, Won J-C, et al. Role of hypothalamic Foxo1 in the regulation of food intake and energy homeostasis. Nat Neurosci. 2006;9:901–6.

    CAS  PubMed  Google Scholar 

  67. 67.

    Kitamura T, Feng Y, Kitamura YI, Chua SC, Xu AW, Barsh GS, et al. Forkhead protein FoxO1 mediates Agrp-dependent effects of leptin on food intake. Nat Med. 2006;12:534–40.

    CAS  PubMed  Google Scholar 

  68. 68.

    Varela L, Horvath TL. Leptin and insulin pathways in POMC and AgRP neurons that modulate energy balance and glucose homeostasis. EMBO Rep. 2012;13:1079–86.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Hahm S, Mizuno TM, Wu TJ, Wisor JP, Priest CA, Kozak CA, et al. Targeted deletion of the Vgf gene indicates that the encoded secretory peptide precursor plays a novel role in the regulation of energy balance. Neuron. 1999;23:537–48.

    CAS  PubMed  Google Scholar 

  70. 70.

    Hahm S, Fekete C, Mizuno TM, Windsor J, Yan H, Boozer CN, et al. VGF is required for obesity induced by diet, gold thioglucose treatment, and agouti and is differentially regulated in pro-opiomelanocortin-and neuropeptide Y-containing arcuate neurons in response to fasting. J Neurosci. 2002;22:6929–38.

    CAS  PubMed  PubMed Central  Google Scholar 

  71. 71.

    Jethwa PH, Warner A, Nilaweera KN, Brameld JM, Keyte JW, Carter WG, et al. VGF-derived peptide, TLQP-21, regulates food intake and body weight in Siberian hamsters. Endocrinology. 2007;148:4044–55.

    CAS  PubMed  Google Scholar 

  72. 72.

    Levi A, Ferri G-L, Watson E, Possenti R, Salton SR. Processing, distribution, and function of VGF, a neuronal and endocrine peptide precursor. Cell Mol Neurobiol. 2004;24:517–33.

    CAS  PubMed  Google Scholar 

  73. 73.

    Grosse J, Heffron H, Burling K, Hossain MA, Habib AM, Rogers GJ, et al. Insulin-like peptide 5 is an orexigenic gastrointestinal hormone. Proc Natl Acad Sci. 2014;111:11133–8.

    CAS  PubMed  Google Scholar 

  74. 74.

    Volkow ND, Wang G-J, Baler RD. Reward, dopamine and the control of food intake: implications for obesity. Trends Cogn Sci. 2011;15:37–46.

    CAS  PubMed  Google Scholar 

  75. 75.

    An H, Ordureau A, Paulo JA, Shoemaker CJ, Denic V, Harper JW. TEX264 is an endoplasmic reticulum-resident ATG8-interacting protein critical for ER remodeling during nutrient stress. Mol Cell. 2019;74:891–908.

    CAS  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Chino H, Hatta T, Natsume T, Mizushima N. Intrinsically disordered protein TEX264 mediates ER-phagy. Mol Cell. 2019;74:909–21.

    CAS  PubMed  Google Scholar 

  77. 77.

    Bouma HR, Carey HV, Kroese FG. Hibernation: the immune system at rest? J Leukoc Biol. 2010;88:619–24.

    CAS  PubMed  Google Scholar 

  78. 78.

    Ebersole JL, Cappelli D. Acute-phase reactants in infections and inflammatory diseases. Periodontol. 2000;23:19–49.

    CAS  Google Scholar 

  79. 79.

    Olivier E, Soury E, Ruminy P, Husson A, Parmentier F, Daveau M, et al. Fetuin-B, a second member of the fetuin family in mammals. Biochem J. 2000;350:589–97.

    CAS  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Ombrellino M, Wang H, Yang H, Zhang M, Vishnubhakat J, Frazier A, et al. Fetuin, a negative acute phase protein, attenuates TNF synthesis and the innate inflammatory response to carrageenan. Shock (Augusta, Ga). 2001;15:181–5.

    CAS  Google Scholar 

  81. 81.

    Delanghe JR, Langlois MR. Hemopexin: a review of biological aspects and the role in laboratory medicine. Clin Chim Acta. 2001;312:13–23.

    CAS  PubMed  Google Scholar 

  82. 82.

    Quaye IK. Haptoglobin, inflammation and disease. Trans R Soc Trop Med Hyg. 2008;102:735–42.

    CAS  PubMed  Google Scholar 

  83. 83.

    Tseng CF, Lin CC, Huang HY, Liu HC, Mao SJ. Antioxidant role of human haptoglobin. Proteomics. 2004;4:2221–8.

    CAS  PubMed  Google Scholar 

  84. 84.

    Stutchfield BM, Antoine DJ, Mackinnon AC, Gow DJ, Bain CC, Hawley CA, et al. CSF1 restores innate immunity after liver injury in mice and serum levels indicate outcomes of patients with acute liver failure. Gastroenterology. 2015;149:1896–909.

    CAS  PubMed  PubMed Central  Google Scholar 

  85. 85.

    Long AT, Kenne E, Jung R, Fuchs TA, Renné T. Contact system revisited: an interface between inflammation, coagulation, and innate immunity. J Thromb Haemost. 2016;14:427–37.

    CAS  PubMed  Google Scholar 

  86. 86.

    Bohr J, Wickbom A, Hegedus A, Nyhlin N, Hörnquist EH, Tysk C. Diagnosis and management of microscopic colitis: current perspectives. Clin Exp Gastroenterol. 2014;7:273.

    PubMed  PubMed Central  Google Scholar 

  87. 87.

    Banoth B, Chatterjee B, Vijayaragavan B, Prasad M, Roy P, Basak S. Stimulus-selective crosstalk via the NF-κB signaling system reinforces innate immune response to alleviate gut infection. Elife. 2015;4:e05648.

    PubMed Central  Google Scholar 

  88. 88.

    Banerjee S, Shah SK, Melnyk SB, Pathak R, Hauer-Jensen M, Pawar SA. Cebpd is essential for gamma-tocotrienol mediated protection against radiation-induced hematopoietic and intestinal injury. Antioxidants. 2018;7:55.

    PubMed Central  Google Scholar 

  89. 89.

    Marakasova ES, Eisenhaber B, Maurer-Stroh S, Eisenhaber F, Baranova A. Prenylation of viral proteins by enzymes of the host: virus-driven rationale for therapy with statins and FT/GGT1 inhibitors. Bioessays. 2017;39:1700014.

    Google Scholar 

  90. 90.

    Zhou Y, Hou Y, Shen J, Huang Y, Martin W, Cheng F. Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2. Cell Discov. 2020;6:1–18.

    PubMed  PubMed Central  Google Scholar 

  91. 91.

    Huang B, Wu P, Bowker-Kinley MM, Harris RA. Regulation of pyruvate dehydrogenase kinase expression by peroxisome proliferator–activated receptor-α ligands, glucocorticoids, and insulin. Diabetes. 2002;51:276–83.

    CAS  PubMed  Google Scholar 

  92. 92.

    Harris R, Bowker-Kinley MM, Huang B, Wu P. Regulation of the activity of the pyruvate dehydrogenase complex. Adv Enzym Regul. 2002;42:249–59.

    CAS  Google Scholar 

  93. 93.

    Wijenayake S, Tessier SN, Storey KB. Regulation of pyruvate dehydrogenase (PDH) in the hibernating ground squirrel, (Ictidomys tridecemlineatus). J Therm Biol. 2017;69:199–205.

    CAS  PubMed  Google Scholar 

  94. 94.

    Andrews MT. Molecular interactions underpinning the phenotype of hibernation in mammals. J Exp Biol. 2019;222:jeb160606.

    PubMed  Google Scholar 

  95. 95.

    Hooper LV, Midtvedt T, Gordon JI. How host-microbial interactions shape the nutrient environment of the mammalian intestine. Annu Rev Nutr. 2002;22:283–307.

    CAS  PubMed  Google Scholar 

  96. 96.

    Donaldson GP, Lee SM, Mazmanian SK. Gut biogeography of the bacterial microbiota. Nat Rev Microbiol. 2016;14:20–32.

    CAS  PubMed  Google Scholar 

  97. 97.

    Carey HV, Duddleston KN. Animal-microbial symbioses in changing environments. J Therm Biol. 2014;44:78–84.

    PubMed  PubMed Central  Google Scholar 

  98. 98.

    Carey HV, Walters WA, Knight R. Seasonal restructuring of the ground squirrel gut microbiota over the annual hibernation cycle. Am J Physiol Regul Integr Comp Physiol. 2013;304:R33–42.

    CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Ma L and Stephen Rossiter for helpful comments on the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (No. 31570378).

Author information

Affiliations

Authors

Contributions

HS and XM conceived the project; HS and JW analysed the data; YP and YX provided tissue samples; HS and XM wrote the draft of the manuscript and YP and XM edited it; XM provided the grant for this project.

Corresponding author

Correspondence to Xiuguang Mao.

Ethics declarations

Ethics approval and consent to participate

Bat sampling procedures were in accordance with the Regulations for the Administration of Laboratory Animals (Decree No. 2 of the State Science and Technology Commission of the People’s Republic of China on November 14, 1988) and approved by the Animal Ethics Committee of East China Normal University (approval number AR2012/03001).

Consent for publication

No applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1.

Principal component analysis (PCA) and clustering of the 30 samples from 15 individuals based on expression data of 14,009 genes showing gene expression distance of all samples.

Additional file 2: Figure S2.

Number of differentially expressed genes (DEGs) identified in each pairwise comparison of the four states by DESeq2 and edgeR, respectively.

Additional file 3: Figure S3

. Results of GO terms and/or KEGG and Reactome pathway enrichment analysis on upregulated DEGs identified in four comparisons. Description of terms or pathways with same color indicates similar function.

Additional file 4: Figure S4.

Results of KEGG and Reactome pathway enrichment analysis on downregulated DEGs identified in five comparisons. Description of terms or pathways with same color indicates similar function.

Additional file 5: Figure S5.

Results of GO terms and/or KEGG and Reactome pathway functional enrichment analysis on genes in TP-M9, TP-M10 and AS-M1 modules identified in WGCNA. Description of terms or pathways with same color indicates similar function.

Additional file 6: Table S1.

Sequencing data collected from RNA-seq. Table S2. Summaries of assemblies generated by TRINITY.

Additional file 7: Table S3.

Details of DEGs identified in each pairwise comparison of four physiological states using edgeR.

Additional file 8: Table S4.

Details of DEGs identified in each pairwise comparison of four physiological states using DESeq2.

Additional file 9: Table S5.

Results of WGCNA module detection of 1614 DEGs and mean log2CPM of DEGs in each of the four states. Table S6. Gene ontology (GO) terms enriched in DEGs identified in each pairwise comparison of states. Table S7. KEGG and Reactome pathways enriched in DEGs identified in each pairwise comparison of states. Table S8. Summary of KEGG and Reactome pathways and GO terms enriched in DEGs in each WGCNA module. Table S9. Details of DEGs included in PPI network of module TP-M9 associated with Torpor. Table S10. Details of DEGs included in PPI network of TP-M10 module associated with Torpor. Table S11. Details of DEGs included in PPI network of module AS-M1 associated with Arousal. Table S12. DEGs upregulated in winter torpor relative to summer active that were shared by the three tissues (intestine, liver and brain) in Rhinolophus ferrumequinum.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sun, H., Wang, J., Xing, Y. et al. Gut transcriptomic changes during hibernation in the greater horseshoe bat (Rhinolophus ferrumequinum). Front Zool 17, 21 (2020). https://doi.org/10.1186/s12983-020-00366-w

Download citation

Keywords

  • RNA-Seq
  • Food intake
  • Immune system
  • Intestine
  • Mammals