Liver Transcriptome Dynamics During Hibernation Are Shaped by a Shifting Balance Between Transcription and RNA Stability

Abstract

Hibernators dramatically lower metabolism to save energy while fasting for months. Prolonged fasting challenges metabolic homeostasis, yet small-bodied hibernators emerge each spring ready to resume all aspects of active life, including immediate reproduction. The liver is the body's metabolic hub, processing and detoxifying macromolecules to provide essential fuels to brain, muscle and other organs throughout the body. Here we quantify changes in liver gene expression across several distinct physiological states of hibernation in 13-lined ground squirrels, using RNA-seq to measure the steady-state transcriptome and GRO-seq to measure transcription for the first time in a hibernator. Our data capture key timepoints in both the seasonal and torpor-arousal cycles of hibernation. Strong positive correlation between transcription and the transcriptome indicates that transcriptional control dominates the known seasonal reprogramming of metabolic gene expression in liver for hibernation. During the torpor-arousal cycle, however, discordance develops between transcription and the steady-state transcriptome by at least two mechanisms: 1) although not transcribed during torpor, some transcripts are unusually stable across the torpor bout; and 2) unexpectedly, on some genes, our data suggest continuing, slow elongation with a failure to terminate transcription across the torpor bout. While the steady-state RNAs corresponding to these read through transcripts did not increase during torpor, they did increase shortly after rewarming despite their simultaneously low transcription. Both of these mechanisms would assure the immediate availability of functional transcripts upon rewarming. Integration of transcriptional, post-transcriptional and RNA stability control mechanisms, all demonstrated in these data, likely initiate a serial gene expression program across the short euthermic period that restores the tissue and prepares the animal for the next bout of torpor.

Keywords: A-to-I RNA editing; ARE-binding proteins; Ictidomys tridecemlineatus; alternative splicing; hepatocyte.

FIGURE 1

Physiological groups analyzed for differential gene expression in liver using RNA-seq and GRO-seq. (A) Schematic shows Tb characteristics of animals in active (SA and SpD) vs. hibernation seasons (IBA, Ent, ET, LT, Ar). RNA-seq data were collected from n = 5 individuals in SA, IBA, Ent, Ar and SpD. Groseq data were collected for n = 3 individuals in SA, IBA, Ent, ET and LT (see Supplementary Table 1 for metadata). (B) Venn diagram shows common and unique sample groups studied by RNA-seq and GRO-seq.

FIGURE 2

RNA-seq reveals changes in the relative abundance of liver steady-state RNAs as a function of hibernation physiology. Two-dimensional scaling plots obtained by unsupervised (A) or supervised (B) clustering of RNA-seq pass-filter reads by random forest; each label represents one individual (groups are as defined in Figure 1). (C) Schematic of seasonal (SpD - > SA - > IBA) and hibernation (IBA - > Ent - > Ar) cycles showing the number of DE genes across each pairwise transition. As indicated in the schematic legend (left), the number of DE genes increased in the state closest to the head or the tail of the arrow are enumerated in black and gray, respectively. Numbers below report the number of DE genes/the number of pass-filter liver genes evaluated. Line plots of (D) 9 liver DE genes exhibiting log2 fold change > 5; (E) Differentially expressed microRNA precursor genes. (F) Boxplot showing the GC-content of DE genes (dots) across the torpor-arousal cycle. Each of the three possible pairwise comparisons comprising the torpor-arousal cycle depicted in panel C are indicated below, genes elevated in each state of the pair are plotted by color (legend). See also Supplementary Data Sheet 3, Supplementary Table 2.

FIGURE 3

Seasonally decreased metabolic genes dominate differential gene expression in hibernator liver. (A) Heatmap summary of RNA-seq co-expression data for 3,120 DE genes in liver. Colors indicate the relative abundance of steady-state RNA for each physiological state. Pattern names (left) and the number of genes (right) in each co-expression cluster are indicated. (B,C) Gene enrichment categories in the DE genes. Two-dimensional scaling plots for gene enrichments based on similarity of terms, using REVIGO (B) and then replotted in panel (C) after segregating by co-expression pattern and adjusting so that circle diameters are proportional to the number of DE genes in each enrichment category. Colors group enrichments to the indicated broad categories: expression = regulation of gene expression, repair = DNA repair; RNA = RNA splicing and processing, and stress = stress response and signaling. In panel (C), co-expression clusters with < 2 significant enrichments are excluded. See also Supplementary Figure 2 and Supplementary Table 3.

FIGURE 4

GRO-seq reveals transcriptional run-on and failure to terminate during torpor. (A) Browser view shows 13-lined ground squirrel (HiC_Itri_2) genomic region with Nxph4 downstream of Lrp1 and upstream of Shmt2, all transcribed from the minus strand, along with GRO-seq coverage from one representative sample from each state across this region and RNAseq from the common SA and Ent samples. Arrows mark paused polymerase near transcriptional start sites (promoter proximal pause) for Lrp1 and Shmt2. (B) Line graphs plot log2 fold change for these three genes in GRO-seq and RNA-seq data. (C) Distribution of fstitch segments called in GRO-seq data for each hibernation state, relative to the position of the 3′ end (at 0, ± 20kb) of annotated genes. The 3′ ends included in the analysis were required to have fstitch annotations called in all 5 states and these annotations were required to overlap the annotated 3′ end. RNA polymerase is known to accumulate at both the 5′ and 3′ ends of transcripts (Core et al., 2008; Glover-Cutter et al., 2008). See also Supplementary Figures 3, 6, and 13.

FIGURE 5

Differentially transcribed genes in liver during hibernation captured by GRO-seq. (A) Two-dimensional scaling plot shows unsupervised random forest clustering of individuals based on the full GRO-seq dataset. (B) Numbers of differentially transcribed genes between pairs of sequential stages in the torpor-arousal cycles of hibernation, compared to summer. As indicated in the legend below, the number of differentially transcribed genes increased in the state closest to the head or the tail of the arrow are enumerated in black and gray, respectively. Numbers below are the number DE/total fstitch gene-body transcription units. (C) Heatmap summarizes the GRO-seq coexpression clusters with the pattern indicated on the left and number of genes in each cluster indicated on the right. See also Supplementary Figure 5 and Supplementary Table 4.

FIGURE 6

Comparison of DE genes identified by RNA-seq and GRO-seq. (A) Venn diagram indicates the number of DE genes in just RNA-seq (blue), just GRO-seq (orange), or both. (B) Common DE genes in pairwise sequential transitions; for this comparison, LT (GRO-seq) and Ar (RNA-seq) are taken as comparable states. (C–F) Scatter plots of log2 fold change in RNA-seq vs. GRO-seq data for all 770 genes in the indicated pairwise transition (see panel B). DE genes for the indicated pair are plotted in dark red, all others are gray. Correlations for the 770 genes across each pairwise transition were: (C) IBA-SA, r = 0.76, p-value < 2.2e-16; (D) IBA-Ent, r = 0.204, p-value = 1.13e-08; (E) Ent-LT/Ar, r = 0.181, p-value = 3.97e-7; (F) IBA-LT/Ar, r = 0.174, p-value = 1.18e-06. (G) Top DE discordant protein coding genes in RNA-seq vs. GRO-seq data.

FIGURE 7

A-to-I RNA editing occurs in the liver during torpor. (A) Summary of single-nucleotide variants detected in RNA-seq data (but not in GRO-seq data) with allele frequencies that are significantly (FDR < 0.05) variable across hibernation stages. (B) Summary of pairwise analysis of A-to-I editing sites with differential editing frequencies between each stage. FDR < 0.05 is considered significant. (C) Euler-diagram comparing the A-to-I RNA editing sites with enhanced editing during torpor in brain tissue (Riemondy et al., 2018) and in liver, p-value < 2.3 × 10^{– 308}. (D) Heatmap of the frequency of editing sites with cold-enriched editing in the liver. Only editing sites with sufficient counts in greater than six samples are shown. (E) Summary of the genomic positions of cold-enriched liver editing sites. (F) Summary of the predicted impacts of A-to-G substitution on mRNA processing and translational activities. See also Supplementary Table 6 and Supplementary Figure 8.

FIGURE 8

Summary of alternative splicing in the liver during hibernation. Heatmaps show summary of (A) mean dPSI, and (B) retained introns relative to SA. Numbers on the right give the number of genes in each of the indicated cluster patterns named on the left. Splice graphs from MAJIQ illustrate temperature dependent alternative splicing for (C) SRSF6 and (D) Trab2b. See also Supplementary Table 7.