Hepatic transcriptome of the freeze-tolerant Cope’s gray treefrog, Dryophytes chrysoscelis: responses to cold acclimation and freezing
BMC Genomics volume 21, Article number: 226 (2020)
Cope’s gray treefrog, Dryophytes chrysoscelis, withstands the physiological challenges of corporeal freezing, partly by accumulating cryoprotective compounds of hepatic origin, including glycerol, urea, and glucose. We hypothesized that expression of genes related to cryoprotectant mobilization and stress tolerance would be differentially regulated in response to cold. Using high-throughput RNA sequencing (RNA-Seq), a hepatic transcriptome was generated for D. chrysoscelis, and gene expression was compared among frogs that were warm-acclimated, cold-acclimated, and frozen.
A total of 159,556 transcripts were generated; 39% showed homology with known transcripts, and 34% of all transcripts were annotated. Gene-level analyses identified 34,936 genes, 85% of which were annotated. Cold acclimation induced differential expression both of genes and non-coding transcripts; freezing induced few additional changes. Transcript-level analysis followed by gene-level aggregation revealed 3582 differentially expressed genes, whereas analysis at the gene level revealed 1324 differentially regulated genes. Approximately 3.6% of differentially expressed sequences were non-coding and of no identifiable homology. Expression of several genes associated with cryoprotectant accumulation was altered during cold acclimation. Of note, glycerol kinase expression decreased with cold exposure, possibly promoting accumulation of glycerol, whereas glucose export was transcriptionally promoted by upregulation of glucose-6-phosphatase and downregulation of genes of various glycolytic enzymes. Several genes related to heat shock protein response, DNA repair, and the ubiquitin proteasome pathway were upregulated in cold and frozen frogs, whereas genes involved in responses to oxidative stress and anoxia, both potential sources of cellular damage during freezing, were downregulated or unchanged.
Our study is the first to report transcriptomic responses to low temperature exposure in a freeze-tolerant vertebrate. The hepatic transcriptome of Dryophytes chrysoscelis is responsive to cold and freezing. Transcriptomic regulation of genes related to particular pathways, such as glycerol biosynthesis, were not all regulated in parallel. The physiological demands associated with cold and freezing, as well as the transcriptomic responses observed in this study, are shared with several organisms that face similar ecophysiological challenges, suggesting common regulatory mechanisms. The role of transcriptional regulation relative to other cellular processes, and of non-coding transcripts as elements of those responses, deserve further study.
Winter in temperate regions poses myriad challenges to vertebrates. While some endotherms endure winter’s low temperatures by elevating metabolic heat production and remaining euthermic, others abandon thermal homeostasis, suppressing metabolism, and tolerating reduced body temperature. For overwintering ectotherms, avoidance of somatic freezing can be achieved by locating protected hibernacula or by supercooling . However, a few temperate-zone ectothermic vertebrates survive winter by tolerating freezing of up to 65% of their body fluids [1, 2]. Somatic freezing generates physiological challenges beyond those encountered by endotherms and non-freezing ectotherms, including osmotic imbalances, compromised oxygen delivery to cells, and the potential for tissue damage from ice. To survive these stresses, freeze-tolerant anurans employ several molecular and physiological strategies [2, 3]. Most freeze-tolerant frogs accumulate cryoprotectants: low molecular weight compounds that reduce ice content, stabilize membranes and proteins, and may serve as metabolic substrates during subzero exposures [2, 4]. Furthermore, although cold acclimation and freezing are associated with metabolic suppression , genes related to production and distribution of cryoprotectants are upregulated during these periods [5,6,7], as are antioxidant enzymes  and the cellular stress response , partially via upregulation of chaperone proteins .
Gray treefrogs, Dryophytes (formerly Hyla) chrysoscelis and D. versicolor , inhabit eastern North America. In fall, treefrogs shift from arboreal microhabitats to sites below terrestrial leaf litter where, in more northern latitudes, they may experience subzero temperatures and somatic freezing [12, 13]. Upon initiation of freezing, treefrogs mobilize cryoprotective glucose and glycerol derived from liver glycogen catabolism [14,15,16,17]. These species also accumulate cryoprotective glycerol and urea during cold acclimation in anticipation of freezing [5, 13, 18,19,20,21] a response similar to that observed in freeze-tolerant insects  and cold-tolerant fish . In D. chrysoscelis, the anticipatory accumulation of cryoprotectant occurs concurrently with changes in kidney function  and adjustments in aquaglyceroporin expression and localization [7, 20], which likely facilitate permeation of these solutes into various tissues. The molecular mechanisms regulating cold-acclimation and freeze-tolerance responses in this species remain largely unknown.
Previous studies have used pathway- and gene-specific approaches to investigate molecular regulation of cold and freezing responses in several organisms [5, 6, 24]. The development of high-throughput RNA sequencing (RNA-Seq) permits a transcriptome-wide investigation of responses to a variety of stimuli, in a way that allows for both testing of hypotheses and discovery of genes and responses integrated across multiple pathways (e.g. [25,26,27]). A transcriptome-level investigation of the hepatic responses of D. chrysoscelis to low temperatures may be particularly informative for understanding cold acclimation and freeze tolerance, as the liver is essential in cryoprotectant mobilization during these challenging periods . In this study, it was hypothesized that hepatic pathways related to cryoprotectant mobilization and stress tolerance would be differentially regulated to meet the physiological challenges encountered during fall and winter. To test these hypotheses, a hepatic transcriptome was generated using RNA-Seq and gene expression was compared among frogs that were warm-acclimated (hereafter called “warm frogs”), cold-acclimated (“cold frogs”), and frozen (“frozen frogs”).
De novo assembly and annotation
A de novo transcriptome was assembled from warm, cold, and frozen frogs (n = 4 for all treatments) to increase transcriptional diversity (Table 1). Twelve cDNA libraries were generated from liver tissue of D. chrysoscelis. These libraries produced a total of 886 million reads, including 159,556 transcripts, and the assembly had an average read depth of 142. Transcripts were on average 676 bp long, ranging from 201 to almost 16,000 bp (Table 1, Fig. 1). FastQC analyses (FastQC, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) confirmed per base sequence quality, per sequence quality, and per tile sequence quality throughout the 150 bp sequence length per read for each biological sample, with zero per base N content reported. Following BLASTX (or TBLASTX as appropriate) searches of the Uniprot Swiss-Prot database  using Trinotate  52,919 transcripts were identified; using BLAST2GO  (BioBam, Spain) an additional BLASTX search of the NCBI vertebrate non-redundant (nr) protein database and of the reference protein sequences (RefSeq) for Xenopus sp identified 62,420 transcripts . Combined, the two approaches identified 62,454 transcripts (Table 2). The taxonomic distribution of the BLAST results revealed that, with Trinotate, 18% of transcripts matched Xenopus sp. and 36% matched Homo sapiens sequences, whereas in Blast2GO, 68 and 1% of the transcripts matched Xenopus sp. and H. sapiens sequences, respectively (Fig. 2). The Japanese treefrog (D. japonicus), a species closely related to D. chrysoscelis, was not represented in the Trinotate results (it was absent from the Swiss-Prot database at the time of search) and in Blast2GO accounted for fewer than 1% of transcripts identified.
A total of 48,865 and 49,989 transcripts were annotated with gene ontology (GO) terms  in Trinotate and Blast2GO respectively, and 53,583 transcripts (34% of the total number of transcripts) were annotated when combining the results from both analyses (Table 2). Trinotate identified 21,065 genes, and 19,705 of these were annotated (Table 2). In contrast, Blast2GO identified 34,936 genes, and 29,879 of these were annotated. To maximize coverage, the transcriptome annotated with Trinotate was used, and the genes that were only identified in Blast2GO were added to that transcriptome. Approximately 0.1% of the total transcripts were identified as non-coding functional RNAs upon searching the Rfam database  using Blast2GO (Table 3). Of note, the transcriptome of D. chrysoscelis included both C/D box and H/ACA box small nucleolar RNAs and 20 different microRNAs: let-7 and miR-7, 10, 15, 16, 22, 23, 24, 27, 29, 30, 101, 122, 142, 144, 145, 194, 214, 221, and 363.
Differential gene expression
Expression patterns of transcripts and genes in warm (n = 4), cold (n = 4) and frozen (n = 4) animals were compared using DESeq2 . We were interested not only in determining which genes showed a significant change (of at least two-fold) in transcriptional output when considering all transcript isoforms, but also which genes had a minimum of one transcript responding to treatment (with at least a two-fold change). We chose this approach as some genes may be identified as differentially expressed without having any of the transcript isoforms show differential expression according to our criteria.
Transcript-level analysis followed by gene-level aggregation revealed 3582 genes that were differentially expressed. In contrast, analysis at the gene level revealed 1324 differentially regulated genes. In both transcript and gene-level analyses 1191 genes were differentially expressed. Using transcript-level analysis, cold acclimation resulted in downregulation of 629 genes and upregulation of 1917 genes, whereas freezing resulted in downregulation of 1093 genes and upregulation of 2223 genes, relative to warm-acclimated animals (Fig. 3). Surprisingly, only 18 genes were upregulated in frozen frogs relative to cold animals, and only 7 were significantly downregulated. Gene-level analysis yielded a lower number of differentially expressed genes in all comparisons (Table 4). When contrasting the two analyses, 86–91% of the genes identified as differentially expressed using the gene-level analysis were also differentially expressed in the transcript-level analysis. However, only 20–44% of the genes identified in the transcript-level analysis were shared by both approaches.
As detailed above, a substantial number of sequences were not annotated due to the absence of homology with known genes or transcripts and contained no predicted functional open reading frame. However, a subset of these non-annotated, non-coding sequences (representing ~ 3.6% of the total transcripts) exhibited differential expression in this study when analyzed with all the transcripts (regardless of annotation status). Within this subset, greater than twice the number of sequences were upregulated as were downregulated in liver from cold and frozen frogs relative to warm (Table 5). This regulatory pattern is consistent with that observed with protein coding transcripts (Table 4), suggesting that both coding and non-coding RNAs are subject to thermal state regulation. From these non-annotated transcripts that were differentially expressed, several non-coding RNA families were identified: several C/D box small nucleolar RNAs (SNORDs), H/ACA box small nucleolar RNAs (SNORAs), and two microRNAs, miR-30 and miR-142, were upregulated during cold acclimation and freezing (Table 6).
Functional classification of transcriptome and differentially regulated genes
Functional classification of the annotated transcriptome regarding Biological Process categorized 44% of genes as related to cellular processes and 38% related to cellular metabolism (Fig. 4). The most common molecular functions were Catalytic (33%) and Binding (27%), and approximately 24% of the genes were associated with the cytoplasm and nucleus (Intracellular) (Fig. 4). We followed a conservative approach and investigated the genes that were differentially expressed based on both gene- and transcript-level analyses (i.e., genes with at least a two-fold change in transcriptional output and with at least one transcript isoform showing a two-fold change). Cold acclimation and freezing resulted in differential expression of genes related to several GO slim Biological Process (BP) categories including nucleotide and nucleic acid metabolism, lipid metabolism, and stress, among others (Fig. 5; Additional file 1: Table S1). Few genes were differentially expressed between cold and frozen frogs, and thus most of the analysis herein will focus on comparisons with warm animals. Genes related to carbohydrate and lipid metabolism were mostly downregulated during cold acclimation and freezing, relative to warm animals (Fig. 5A and B). Among these were glyceraldehyde-3-phosphate dehydrogenase (GAPDH), glycerol kinase (GLPK), acyl-CoA dehydrogenase family member 10 (ACAD10) and diacylglycerol O-acyltransferase 2 (DGAT 2) (Additional file 1: Table S1). However, some genes involved in these same general processes, such as hexokinase (HK2) and pyruvate dehydrogenase (PDK2), were upregulated (Additional file 1: Table S1). Most genes associated with nucleotide and protein metabolism were upregulated, including proteins associated with ubiquitination, although some genes involved in nucleotide and nucleic acid metabolism, like elongation factor 1-alpha 1 (EEF1A1) and signal transducer and activator of transcription 2 (STAT2), as well as genes involved in amino acid metabolism, including alanine aminotransferase 2 (GPT2) and glutamate dehydrogenase (GLUD1), were downregulated in cold and frozen animals (Additional file 1: Table S1). Expression of several transport-related genes was generally higher in cold and frozen than in warm animals, but a few notable exceptions included the solute carrier family 2, facilitated glucose transporter members 8 and 10 (SLC2A8, SLC2A10), and sodium/glucose cotransporter 4 (SLC5A9) (Additional file 1: Table S1), all of which were downregulated in the cold. Stress-related genes were differentially regulated under cold acclimation (Fig. 5a and b). Specifically, heat shock proteins 90 (HSP90AA1) and 70 (HSP70) were upregulated in cold frogs, whereas catalase (CAT) and caspase 3 (CASP3) were downregulated in cold and frozen animals, relative to the warm group (Additional file 1: Table S1). Comparing frozen and cold frogs directly revealed only 3 genes with differential expression: complement component 1Q subcomponent-binding protein, mitochondrial; F-box/LRR-repeat protein 5; and the gene for veficolin-1 (Additional file 1: Table S2). Among genes that were upregulated in cold animals, there was an approximate two-fold significant enrichment (FDR < 0.05) of genes associated with RNA metabolism and with nitrogen compound metabolism. On the other hand, genes that were downregulated in cold frogs were enriched for processes associated with organic acid metabolism and amino acid metabolism, a pattern also detected in frozen frogs, relative to warm animals (Fig. 6a and b).
Temperature is one of the most pervasive factors influencing the physiology of ectotherms. Cope’s gray treefrog, Dryophytes chrysoscelis, undergoes molecular, cellular, and physiological adjustments during cold exposure that ensure functional integrity during corporeal freezing . In addition to the metabolic suppression induced by low temperatures, physiological challenges associated with internal ice formation include dehydration and anoxia, along with possibilities of oxidative stress and physical damage from ice [2, 3]. The liver plays a paramount role in anuran cryoprotectant mobilization, a process essential for freezing survival [2, 3]. The molecular processes involved in hepatic responses to low temperature and freezing are complex and not completely understood. The present study investigated responses of the hepatic transcriptome to cold acclimation and freezing in Dryophytes chrysoscelis, finding that several cellular and metabolic processes are transcriptionally responsive to thermal changes, potentially promoting cryoprotectant mobilization, stress responses, and transcriptional repression. The environmental and organismal changes that occur during the shift from warm to cold condition— low temperature, anoxia, dehydration, reduced photophase, and aphagia— may have promoted these responses [2, 3].
The hepatic transcriptome of D. chrysoscelis
The RNA-Seq analysis undertaken in this study relied on a de novo transcriptome assembly, generating a high number of transcripts and an N50 on par with those in other RNA-Seq studies that involved a de novo assembly in amphibians (anurans or urodeles) [27, 35, 36]. Two different approaches were used for transcriptome annotation: one approach used Trinotate, which relies on Swiss-Prot as the main database for BLASTX searches; the other relied on Blast2GO to BLASTX sequences on NCBI’s nr protein database, as well as Xenopus sp. RefSeq. Blast2GO identified more transcripts than Trinotate, and the taxonomic distribution of the BLAST hits differed between the two approaches, reflecting the inherent bias of species’ representation in the Swiss-Prot database (which contains more than 20,200 H. sapiens sequences and approximately 6300 anuran sequences, but the latter representing few genera other than Xenopus).
Approximately 34% of transcripts in the present assembly were annotated, within the range of values (32–57%) reported in similar studies [27, 35,36,37]. In D. chrysoscelis, the most common Metabolic Function category in the transcriptome was Catalytic, as also was observed in the transcriptome of the green-striped burrowing frog (Cyclorana alboguttata) . In contrast, in the crab-eating frog (Fejervarya cancrivora) and the green odorous frog (Odorrana margaretae), Binding was the most common metabolic function [27, 35]. As in other studies reporting amphibian transcriptomes [27, 35, 36], Cellular and Metabolic Processes were the most abundant Biological Process categories in the hepatic transcriptome of D. chrysoscelis. The present study focused on the hepatic transcriptome, whereas other studies have analyzed and reported a merged gene ontology classification of multiple tissues [27, 35, 36]; despite the distinct approaches and evolutionary histories of the species analyzed, the overall composition of the transcriptomes of D. chrysoscelis and these anurans is similar.
In addition to protein-coding genes, we also identified transcripts representing several categories of non-coding RNAs. These included several small nucleolar RNAs (snoRNA); 20 different microRNAs, several of which (miR-16, 24, 30, 145, 145, 214, and 363) were also identified in other freeze tolerant amphibians [38, 39] and a number of putative long non-coding RNAs.
Differential gene expression in liver of D. chrysoscelis
Analysis of differential transcript and gene expression was conducted on the 34% of the transcriptome that was annotated. Our criterion for differential gene expression was fairly conservative, as we did not only want to detect genes that had at least one transcript showing differential expression, but that also had an overall change in transcriptional output . In response to cold acclimation and freezing, 5.1% of the annotated genes (5.6% under the gene-level analysis and 15.3% under the transcript-level analysis) were differentially expressed in D. chrysoscelis. Collapsing alternative transcripts into a gene-based analysis allowed for a more direct comparison with other published studies. However, in doing so, the gene-level analysis underrepresents the complexity of acclimation-state dependent differential gene expression. Transcript-level analyses revealed 2–4 fold the number of differentially expressed sequences relative to those identified through the gene level analysis (Table 4). The identification of alternative transcripts, and therefore, potentially different protein isoforms that are differentially expressed suggests a role for alternative mRNA processing as a component of acclimation-state dependent gene regulation.
Most differentially expressed genes were found when comparing either cold or frozen frogs to warm animals. Interestingly, 56% of differentially expressed genes (75% of differentially expressed transcripts) showed an increase in expression in cold- vs. warm-acclimated animals. The number of differentially expressed genes and transcripts was further increased in frozen animals relative to warm-acclimated treefrogs, though the number of up- and downregulated genes was similar; in contrast the number of upregulated transcripts was double that of those downregulated, relative to warm-acclimated animals. Only a few genes were differentially expressed between cold and frozen frogs, regardless of the analysis used. These results suggest that physiological adjustments required to survive freezing are initiated during cold acclimation, but that, under conditions of low temperature and metabolic suppression at the time of ice formation, further genomic response is limited. Although a longer duration of freezing might have induced additional transcriptomic adjustments, transcription and translation of new proteins may proceed slowly enough at these low temperatures that other, less energetic mechanisms, such as epigenetic modifications, signal transduction, and actions of microRNAs, predominate . Due to the few differences in gene expression between cold and frozen frogs, the remaining discussion focuses mostly on the comparisons of cold and frozen frogs vs. warm animals.
Differential expression of genes involved in cryoprotectant mobilization
Glycerol is one of the major cryoprotectants in D. chrysoscelis, and this metabolite can be accumulated in the fall , and upon initiation of freezing , likely as a result of hepatic glycogen catabolism [15, 17]. In this study we found evidence of transcriptional regulation of glycerol mobilization pathways, as enzymes associated with carbohydrate metabolism, and glycerol metabolism in particular, showed differential expression, a pattern found at the transcript and gene-level analyses. Counter to our expectations, that response did not entail enhanced expression of genes promoting glycerol synthesis. Those enzymes, including glycerol-3-phosphate dehydrogenase and glycerol-3 phosphatase, did not show differential gene expression under cold or frozen conditions. Further, in the pathway for glycerol biosynthesis, transcript levels of the key enzyme glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were diminished in cold and frozen animals. Downregulation of this gene, along with decreased gene expression of glutamate dehydrogenase (GLUD1) and alanine aminotransferase (GPT2), suggests a reduced role for glyceroneogenesis in glycerol accumulation in cold and frozen frogs . At the same time, expression of glycerol kinase (GK), an enzyme essential for the phosphorylation and metabolism of glycerol in liver, was decreased in cold and frozen animals. Downregulation of GK transcription likely decreases glycerol metabolism and thereby contributes to its accumulation for cryoprotective purposes. Curiously, expression of another member of the glycerol kinase family, glycerol kinase 5, was elevated in cold and frozen animals. The function of this kinase in treefrogs is currently unknown but in mice the enzyme is restricted to the skin and has some glycerol kinase activity . A decrease in expression of genes associated with amino acid metabolism is further supported by the GO enrichment analysis which shows downregulated genes being enriched in amino acid metabolism processes. Collectively, these findings suggest that glycerol accumulation during cold acclimation and freezing does not involve transcriptional upregulation of hepatic glycerogenesis, but rather a combination of other regulatory mechanisms, potentially including diminished glycerol metabolism, provision of glycerol from non-hepatic sources like fat bodies, or regulation of transport and distribution of the metabolite.
Transmembrane transport of glycerol occurs largely via aquaglyceroporins, a subclass of the aquaporin protein family that confers permeability to glycerol and water [5, 20]. In previous studies of D. chrysoscelis, cold acclimation induced increases in hepatic transcript levels of HC-1 (a water-permeable aquaporin) and HC-3 (an aquaglyceroporin permeable to water and glycerol), and these increases accompanied a cold-induced increase of glycerol in plasma . Hepatic protein levels of HC-9, another glycerol-permeable aquaglyceroporin, also increased during freezing in D. chrysoscelis , without an associated change in transcript levels. We did not detect changes in transcript levels of either aquaglyceroporin (HC-3 or HC-9) in the current study. The contribution of aquaglyceroporins to glycerol mobilization (export from hepatocytes, uptake in other tissues) may occur via post-transcriptional mechanisms, such as trafficking of the proteins between cytoplasm and plasma membrane.
In treefrogs, mobilization of cryoprotective glucose takes place when freezing begins [14, 15], a process likely triggered by an enzymatic cascade that results in glycogen breakdown and glucose accumulation . Glycogen breakdown in D. chrysoscelis is transcriptionally promoted by an increase in expression of phosphorylase kinase b (PHKA2, PHKG2) in frozen frogs, and of the β catalytic subunit of protein kinase A (PRKACB), which may further potentiate glycogenolysis and subsequent accumulation of glucose and glycerol [4, 6]. In D. chrysoscelis, glucose export may also be promoted transcriptionally through the upregulation of glucose-6-phosphatase (G6PC), which possibly increases glucose available for hepatic export . Glucose accumulation and export is potentially further augmented by the increased expression of pyruvate dehydrogenase kinase 2 (PDK2), which regulates the pyruvate dehydrogenase complex thus inhibiting the oxidative decarboxylation of pyruvate to acetyl-CoA, as well as downregulation of glycolytic enzyme genes like PFKM and GAPDH. The potential downregulation of glycolysis at the level of the hepatic transcriptome reinforces the enzymatic suppression of this process, a response previously reported in other anurans . Curiously, expression of glucokinase regulatory protein gene (GCKR) decreased in cold and frozen frogs, and if reflected at the protein level, would favor glycolysis, a pattern that may be reinforced by the increase in hexokinase (HK2) also detected in cold and frozen frogs.
Facilitative glucose transporters (GLUTs), members of the SLC2A gene family, have an important role in ensuring the export and import of glucose in various tissues and possibly contribute to the flux of hepatic cryoprotective glucose during freezing. Transcript levels of GLUT 8 and GLUT 10 decreased during freezing in D. chrysocelis; however, the role of these transporters in anuran glucose homeostasis is unknown. In contrast, GLUT2, a major hepatic glucose transporter in vertebrates, was not responsive to low temperature exposure or freezing in this study, although protein and transcript levels of GLUT2 increase during freezing in R. sylvatica . The absence of a freezing response by GLUT2 transcripts in D. chrysoscelis may be in part responsible for the low hepatic export of glucose we detected in a previous study . In addition, SGLT4, a Na+-dependent glucose transporter (member of the SLC5A gene family) was downregulated in frozen animals, compared with the warm condition, suggesting further reduction of hepatic import of glucose during a period when glucose export is being maximized.
Overall, our results suggest a complex, multi-level control of expression of transcripts governing cryoprotectant accumulation, involving regulation of genes involved in cryoprotectant synthesis, metabolism, and transport. The temperature-sensitivity of these transcriptional processes, and also of those involved in synthesizing and processing the proteins that derive from those transcripts, ultimately determines the seasonal pattern of cryoprotectant accumulation. The timing of regulated changes in gene expression may well reflect the ability of various biochemical processes—transcription, translation, transport, enzyme-mediated metabolism, and others—to respond adequately to physiological demands at temperatures near or below freezing.
Differential expression of genes involved in stress response
Low temperature exposure and freezing result in drastic changes in the physiology of organisms, triggering a series of reactions common to the cellular stress response . In this study, several pathways related to cellular damage showed differential regulation under cold and freezing exposure. Heat shock proteins 70 (HSP70), 90 (HSP90AA1), 105 (HSPH1), as well as the 78 kDa glucose-regulated protein HSPA5, were upregulated in cold and frozen frogs. Chaperone proteins function in translocation, folding and assembly, and targeting of misfolded proteins for degradation . Although the HSP response has not been previously investigated during freezing in freeze-tolerant anurans, physiological challenges that occur during freezing, including low temperature, anoxia, and dehydration, do trigger upregulation of HSPs in several species of ectotherms [44, 45], including other freeze-tolerant frogs .
The cessation of circulation and ventilation during freezing eventually renders the hepatic tissue anoxic . However, in this study, levels of hypoxia-inducible factor 1 subunit alpha (HIF1A), which plays an essential role in hypoxia response  and is upregulated in heart of frozen R. sylvatica , did not change during freezing (data not shown), perhaps reflecting tissue-specific responses.
Low temperature exposure and freezing are likely to trigger oxidative stress and formation of reactive oxygen species, not only due to metabolic suppression, but also during oscillations in tissue oxygen levels resulting from freezing and thawing . Surprisingly, we detected decreases in expression of genes encoding catalase (CAT), glutathione S-transferase (GSTO1), microsomal glutathione S-transferase 1 (MGST1), and nuclear factor erythroid 2-related factor 2 (NFE2L2), all proteins that promote responses to oxidative stress , suggesting a potential for reduced oxidative stress defenses during cold acclimation and freezing. Furthermore, glutathione-specific gamma-glutamylcyclotransferase 1 (CHAC1), which has been shown to increase oxidative stress by degrading glutathione, was upregulated in cold and frozen frogs. Several organisms upregulate antioxidant responses seasonally or in response to cold exposure [48, 49]. Antioxidant enzyme levels have not been investigated in D. chrysoscelis but may be constitutively elevated as reported in other freeze-tolerant frogs .
Cold exposure and freezing also modulated the hepatic transcriptome of D. chrysoscelis by promoting anti-apoptotic and DNA-repair responses. Caspase 3 (CASP3) which is involved in apoptosis , was downregulated in cold and frozen frogs, whereas the dual specificity protein phosphatase 10 gene (DUSP10), which inhibits p38 MAP kinase and contributes to regulation of apoptosis , was also upregulated in both cold and frozen animals. GADD45, a protein that promotes growth arrest and DNA repair, is responsive to hypoxia , hyperosmotic stress , and is also involved in apoptosis , was upregulated in cold and frozen frogs.
Differential expression of genes related to nucleotide and nucleic acid metabolism
Regulators of transcription and translation were responsive to cold and freezing. Transcription factors can influence transcription by interacting with other transcription factors and RNA polymerase, and by changing the structure of chromatin . The downregulation of several transcription factors, including SOX30 and NFE2L2, suggests a possible mechanism for transcriptional repression. Furthermore, downregulation of the transcription activator STAT2 , and EEF1A1, a key participant in translation and regulator of proteolysis, indicates other potential mechanisms for regulating transcription and translation beyond changes in transcription factors. Decreased levels of EEF1A1 protein have also been observed in cold-exposed X. laevis  and in winter-acclimatized Rana sylvatica .
Differential expression of genes involved in protein, amino acid, and lipid metabolism
Levels of several transcripts involved in the ubiquitin proteasome pathway, including PJA2, RNFL30, UBE2D1, and RBBP6 were upregulated in cold and frozen frogs, whereas FBX15 expression also increased from the cold to the frozen states. The conjugation of ubiquitin with target proteins promotes binding to the proteasome and subsequent protein degradation. Protein ubiquitination is transcriptionally enhanced in cold-exposed fish [48, 56, 57], and frogs , possibly due to cold-induced denaturation , and may serve as a marker for cold exposure . Protein turnover in cold D. chrysoscelis could supply the cell with amino acids to be used for synthesis of ATP, carbohydrates, and proteins, or these might serve as cryoprotectants  or as sources of cryoprotective urea . Given the low oxygen availability in frozen frogs , ATP synthesis from sources other than carbohydrates is unlikely in this state. Moreover, transcripts involved in amino acid trafficking (e.g. GPT2 and GLUD1), and carbohydrate flux (GAPDH) were downregulated in cold and frozen animals, suggesting amino acids derived from proteolysis may be accumulating in the liver, warranting further investigation of hepatic protein metabolism during this period.
Several transcripts associated with β-oxidation of fatty acids were downregulated in our cold and frozen animals, including ACAD10, EHHADH, and DECR1, as well as transcripts of DGAT1, which is essential for triglyceride synthesis. Downregulation of these transcripts may be related to the overall metabolic depression occurring in cold animals, as it curtails unnecessary energetic expenditures incurred by transcription . Furthermore, the predominance of anaerobic metabolism during freezing precludes lipid use for ATP synthesis . Surprisingly, expression of ACACB, involved in fat storage, and DGAT2 increased in cold and frozen frogs.
Differential expression of non-coding RNAs
Approximately 3.6% of the differentially expressed sequences identified in this study were non-coding and showed no homology to sequences in existing databases. In the absence of a reference genome for sequence alignment, it is not possible to determine how many unique transcripts are represented by the 5800 sequences that exhibited differential regulation. The effect of cold exposure and freezing on levels of non-coding RNAs has been investigated in multiple taxa, with particular focus on certain miRNAs [38, 39]. Cold exposure and freezing increased the expression of miR-30 and miR-142 in D. chrysoscelis. In frozen R. sylvatica, miR-30 levels increased with freezing in skeletal muscle ; this microRNA has several demonstrated effects, including inhibiting apoptosis  and protecting against liver fibrosis , but its role in freeze tolerance remains to be elucidated. In contrast, most miRNAs that were differentially regulated in the brain of R. sylvatica had reduced expression in frozen animals, which may suggest a neuroprotective role for those molecules .
Several snoRNAs were responsive to freezing (Table 6). SnoRNAs may have important roles in regulating transcription and translation by modifying rRNA, and by functioning like miRNA , but their role in cold acclimation and freezing is unknown. SnoRA5, which was upregulated in frozen frogs, is associated with kidney reperfusion injury , and its upregulation in frozen animals may be related to the cessation of circulation that occurs in that condition .
In addition to the identified non-coding RNAs, it is likely that among the transcripts in this subset are long non-coding RNAs (lncRNAs) that may contribute to the process of freeze tolerance through mechanisms that remain to be defined. LncRNAs are broadly characterized as non-conserved, transcribed RNAs > 200 bp in length that can reside in the nucleus or in the cytoplasm and can interact with chromatin as well as ribosomes. LncRNAs regulate gene expression through a diversity of mechanisms, including gene silencing, epigenetics, transcriptional and post-transcriptional regulation, RNA processing and transport (reviewed in ). Moreover, tissue and condition specific differential expression of lncRNAs has been documented in a number of biological and pathophysiological conditions, indicating that regulating the expression of lncRNAs is important to function . Sequence analysis of the North American bullfrog genome identified 6223 putative lncRNAs, a subset of which exhibited thyroid responsive DGE, indicating that lncRNAs are present across the anuran order and are subject to differential regulation . LncRNAs have been implicated in cold and freeze response mechanisms in both plant and animal systems. Cold-dependent differential expression and differential alternative splicing of lncRNAs have been shown to correlate with acclimation and freeze tolerance in Arabidopsis [66, 67]. Moreover, differential expression of lncRNA transcripts and associated target genes involved in cold response and apoptosis in frozen and thawed giant panda (Ailuropoda melanoleuca) sperm exhibited co-expression regulatory patterns, suggesting a role for the interaction between lncRNAs and their target mRNAs in cryopreservation, and as an extension, freeze tolerance . Experiments designed to determine the cis- and trans-targets of these putative lncRNAs, and to characterize the potential differentially expressed genes and co-regulation patterns of lncRNAs and target transcripts involved, are planned.
Conservation and specialization of transcriptomic regulation in the cold
Transcriptomic analyses have been used to identify gene regulatory pathways involved in cold and freeze tolerance in other organisms, including hibernating mammals, plants, fishes, insects and other invertebrates. From these studies, diverse functional pathways and processes subject to cold-responsive control have been identified. Hibernation, characterized by alternating periods of euthermia and hypothermia in endotherms, appears to have a complex of both common and species-specific responses. For example, transcriptomic studies in skeletal muscle, adipose tissue and liver from hibernating grizzly bears (Ursus arctos horribilis) revealed differential expression for genes involved in aspects of metabolism (decrease in insulin responsive pathways, upregulation of muscle anabolic pathways), as well as changes in the expression of lncRNAs, in addition to tissue-specific regulatory responses . Thirteen-lined ground squirrels (Ictidomys tridecemlineatus) also differentially express genes in skeletal muscle during hibernation ; however, the particular pathways that are regulated (e.g., protein turnover pathways) differ from those in grizzly bear. Likewise, during discrete phases of their hibernation cycle, white adipose tissue from dwarf lemurs (Cheirogaleus medius) differentially expressed genes related to metabolism, feeding, circadian rhythms and coagulation , a suite of responses that differs from that in adipose tissue of grizzly bears. These findings suggest hibernators may have evolved both common as well as tissue- and species-specific strategies for survival, perhaps related to aspects of body size, metabolic strategy, and natural and phylogenetic history.
In ectotherms, too, transcriptional regulation during cold conditions includes both responses that are common across animals and those that are tissue- or species specific. In threespine stickleback (Gasterosteus aculeatus), cold acclimation resulted in an upregulation in skeletal muscle genes involved in RNA and protein processing and cell cycle regulation and downregulation of transcripts encoding extracellular matrix, angiogenic, and contractile proteins . Analyses of gill and liver transcriptomes from cold-tolerant and cold-sensitive blue tilapia (Oreochromis aureus) showed similar patterns of down regulation of cell-matrix interactions and upregulation of proteolytic enzymes, as well as opposing patterns of upregulation of glycolysis/gluconeogenesis in the gill and amino acid biosynthesis in the liver of cold-sensitive fish, and downregulation of these processes in cold-tolerant fish . More diversity of response is apparent if one extends the comparison to invertebrates (e.g., Drosophila melanogaster, in which cold acclimation induced differential regulation in one third of the transcriptome and one half of the metabolome ) or plants (in which pathways for biosynthesis of secondary metabolites, hormone signal transduction, plant-pathogen interaction, and stress responses all are involved in cold-tolerance [75,76,77]).
In the treefrogs we studied, we identified few differences in gene expression between cold and frozen animals. Nevertheless, there may be categorically important differences between these states in other circumstances or organisms. For example, freeze tolerance in the terrestrial worm (Enchytraeus albidus), is characterized by differential expression of genes involved in anion transport, metabolism, oxidative and metabolic stress reduction, cryoprotectant solute transport and signal transduction . Similarly, differentially expressed genes involved in freeze tolerance of the cricket Gryllus veletis are associated with metabolism, cryoprotectant solute transport, cytoskeleton stability, oxidative and metabolic stress reduction, signal transduction, and cellular processes . It is apparent that the detailed responses to cold and freezing differ among organisms, or even among strains of organisms, and also reflect specialized or condition-specific gene regulatory cascades that are unique and specialized [78, 80, 81]. As described above, transcripts for several genes related to heat shock protein response, DNA repair, and the ubiquitin proteasome pathway showed an increase in abundance in cold and frozen frogs, whereas transcripts for genes involved in oxidative stress and anoxia response, both potential sources of cellular damage during freezing, showed either a reduced abundance or remained unchanged. These transcriptomic results are consistent with correlative trends identified in hepatoproteome studies of thermal responses in R. sylvatica . Interestingly, similar proteomic responses were observed for a subset of genes related to glycogen synthesis and carbohydrate metabolism in response to low temperature exposure (5 °C) in the non-freeze tolerant anuran Xenopus laevis . Beyond transcriptional control, it is also clear from other studies in R. sylvatica that thermal condition dependent regulation of protein expression is influenced by tissue-specific regulation and post-translational modification [82, 83]. Further studies will be required to more broadly define the relation between transcriptomic regulation and protein expression as animals freeze and thaw.
Dryophytes chrysoscelis survives winter’s cold by tolerating freezing of its body fluids, a strategy that requires both a preparatory cold acclimation and a series of cryoprotective responses triggered upon ice nucleation. Hepatic cryoprotectant production during cold acclimation and freezing is of paramount importance to ensure freezing survival. This was the first study to employ RNA-Seq to analyze the freezing response in a freeze-tolerant vertebrate. This study demonstrates that the hepatic transcriptome of Dryophytes chrysoscelis is responsive to cold and freezing with differential expression of numerous genes that contribute to cryoprotectant mobilization, stress response, and transcriptional regulation. The multiple physiological challenges associated with freezing exposure, as well as the transcriptomic responses observed in this study, are shared with several organisms that face similar ecophysiological circumstances [3, 23, 25, 44], suggesting a common regulatory framework. The extent to which these changes in mRNA expression are reflected at the level of the proteome, and the contributions of those changes relative to other cellular responses, remain to be defined. Intriguingly, this study also identified numerous non-coding RNAs that were differentially expressed in response to cold and freezing; the roles of those molecules deserve further study.
Twelve male D. chrysoscelis were collected from ponds in southwest Ohio, USA, within a ten mile radius of Caesar Creek State Park Wildlife Area (39.5065512,-84.0604201) during summer months (May–July). Wild caught frogs (age unknown, weight ~ 7 g) were collected under a permit from the Ohio Division of Wildlife. Frogs were transported to the laboratory at Wright State University and were kept in plastic containers with free access to water. Initially, frogs were warm-acclimated to 22 °C under a 12:12 h light-dark regime, fed crickets thrice weekly, and provided water ad libitum. A random group of these frogs (n = 4) was sampled after a minimum of 8 weeks in acclimation (herein called “warm frogs”). In fall, a random subset of animals was transferred to a cold room and cold-acclimated by being progressively cooled to 5 °C, with an accompanying shift to a 8:16 h light-dark regime, over a 2-month period as previously described . During cold acclimation, each frog was fed until it no longer accepted food, which occurred once the acclimation temperature was 8 °C and lower. In winter, a random subset of the cold-acclimated frogs was sampled (n = 4) (hereafter “cold frogs”) while another subset was used in freezing experiments. If a frog exhibited characteristics outside the normal observed characteristics, the frog was removed from the study in consultation with a laboratory animal research veterinarian and given treatment. Characteristics that would be cause for removal include atypical body posturing, outstretched legs, poor response to stimulus (depending on acclimation state) or other atypical physical characteristics that could be considered symptomatic of disease. Animals to be frozen were placed in individual plastic containers in an incubator set at 2 °C. The temperature of the incubator was then lowered to − 2 °C at a rate of 1 °C day− 1. Once the incubator reached − 2 °C, the unfrozen frogs were individually placed on top of moist gauze and ice was added to the individual containers to ensure crystallization. Immediately after ice was added to the containers, the incubator temperature was lowered to − 2.5 °C. Frogs were visually inspected 24 h post-inoculation for signs of internal ice formation (stiff body and limbs, presence of subcutaneous ice crystals noticeable in the ventral and dorsal surfaces) and sampled shortly thereafter (n = 4) (“frozen frogs”). Internal ice formation was confirmed during tissue sampling by the presence of large pieces of ice surrounding digestive organs, and sheets of ice between the skin and skeletal muscles.
Liver sampling and total RNA isolation
Animals were euthanized following approved WSU Laboratory Animal Care and Use Committee (LACUC) protocols. Warm-acclimated animals were anesthetized by being placed in a shallow volume of skin permeable 0.1% Tricaine-S (MS-222; Pentair Aquatic Systems, Cary NC) in phosphate buffered saline, a well-documented and routinely used anesthetic for amphibians. The anesthetic was taken up by the frog through the pelvic patch. Animals were monitored as they became lethargic and lowered their body further. They were then decapitated before pithing the spinal column. Cold-acclimated animals and frozen animals were directly decapitated and pithed. Warm-acclimated animals were euthanized and dissected at room temperature (22 °C) in an LACUC-approved laboratory space. Cold-acclimated animals and frozen animals were euthanized and dissected in a 4 °C cold room. Animals were euthanized between 10:00 am and 4:00 pm, local EST. The liver was excised and immediately frozen in liquid N2. Samples were stored at − 80 °C until RNA isolation. Total RNA was isolated from 10 to 20 mg of liver using the Qiagen RNAeasy Mini Kit (Qiagen, Hilden, Germany). Tissue samples were homogenized with a hand-held homogenizer with a rotating blade, in isolation buffer, and RNA was isolated following the manufacturer’s instructions. Extracted RNA was treated with DNase I (Ambion, Thermo Fisher Scientific, Waltham, MA, USA) to eliminate potential genomic DNA contamination and precipitated using lithium chloride (Sigma, St. Louis, MO, USA). RNA concentration and quality were determined with a nanophotometer (Denville Scientific Inc., Holliston, MA, USA). Total RNA from 12 male frogs (n = 4 per treatment) was submitted to Cofactor Genomics (http://cofactorgenomics.com/, St. Louis, MO) to be used in high-throughput RNA-Sequencing. Sample size determination was based on statistical analyses of prior gene expression results . As the frogs were wild caught and their ages and degree of inter-relatedness were unknown, it was anticipated that inter-individual variability, independent of acclimation state, but related to age, acute physiological condition, and/or genetic variability, would be observed. Thus, a sample size that slightly exceeds the power analysis generated from prior results was used in this experimental design to address this predicted inter-individual variation.
mRNA purification and next generation sequencing library preparation
Preparation of mRNA library and Next Generation Sequencing were performed by Cofactor Genomics. Samples used for library preparation had an RNA integrity number (RIN) > 7, a UV absorption ratio A260/A280 > 1.9, and A260/A230 > 2. Total RNA samples were incubated with mRNA capture beads to remove contaminating ribosomal RNA using the Kapa Stranded mRNA-Seq kit (Kapa Biosystems, Roche, Pleasanton, CA) following manufacturer’s instructions. The resulting poly(A)-captured mRNA was then fragmented. First-strand cDNA synthesis was performed using random primers and reverse transcriptase in the presence of Actinomycin D, followed by second-strand cDNA synthesis with Rnase H and DNA polymerase I. Double-stranded cDNA was end-repaired and A-tailed for subsequent adaptor ligation. Indexed adaptors were ligated to the A-tailed cDNA. Enrichment by PCR was performed to generate the final cDNA sequencing library. Libraries were sequenced as single-end 75 base pair reads on an Illumina NextSeq500 following the manufacturer’s protocols.
Transcriptome assembly and annotation
Transcriptome assembly and annotation were performed by Cofactor Genomics. Raw sequence data (in FASTQ format) were assessed for quality using FastQC and ribosomal RNA content (sortmeRNA, http://bioinfo.lifl.fr/RNA/sortmerna/). A transcriptome reference was created through the de novo assembly of adapter- and quality-filtered reads. As a foundation for the assembly, Cofactor’s pipeline used the Trinity Suite (Broad Institute, version r20140717) with optimized parameters, resulting in transcript models for alternatively spliced isoforms . For each transcript or patch, the number of mapped reads was determined using Kallisto  and used for statistical analysis.
The de novo assembly was annotated by searching the Uniprot Swiss-Prot database via BLASTX (or TBLASTX when appropriate), following the Trinotate pipeline, with an expectation value (e-value) of 1 × 10− 5 . The transcriptome was also annotated using Blast2GO (BioBam, Spain). Transcripts were imported into Blast2GO and NCBI BLASTX was used to search each transcript against the NCBI vertebrate, non-redundant (nr) protein database with an e-value of 1 × 10− 3, and against the reference protein sequences (RefSeq) for Xenopus sp. The highest scoring BLAST hit was used to assign a gene ID to each transcript. The number of genes identified was estimated by determining the number of distinct gene IDs (Uniprot IDs and NCBI accession numbers) in the transcriptome; transcripts assigned to distinct identifiers were considered different genes. For transcripts with a BLAST match, GO terms were retrieved using the GO-mapping tool (e-value of 1 × 10− 6) and InterProScan tool in Blast2GO, and the annotations from these two tools were merged. To ensure maximal transcriptome coverage, transcriptome annotation performed with Trinotate was merged with the annotation performed in Blast2GO. The remaining unidentified transcripts were searched in the Rfam database, which contains a collection of RNA sequence families of structural RNAs such as non-coding RNA genes, using Blast2GO.
Differential expression analysis
The abundance of annotated transcripts generated in Kallisto were uploaded to R using the package tximport . Differential gene expression was determined at the gene level, and at the transcript level by selecting gene-level summary of the Kallisto data . Differentially expressed genes were also identified from transcript-level data by evaluating if at least one transcript isoform of a gene was differentially expressed, followed by aggregation of transcripts at the gene level . Differential expression among samples was determined using the DESeq2  package for R  in Rstudio (v. 1.0.136). Data were normalized with DESeq2 algorithms. Transcripts with less than one raw count were excluded, and the default DESeq2 algorithms were used to remove outlier transcripts based on Cook’s distances . The Benjamini and Hochberg algorithm  was used to control the false discovery rate (FDR) and transcripts were considered differentially expressed if the expression between two treatments differed by at least |2|-fold with FDR < 0.05.
Differential transcript expression was also determined for the entire transcriptome (including non-annotated transcripts) to determine if non-coding RNAs showed differential expression. The abundance of all transcripts generated in Kallisto was uploaded to R using the package tximport and the analysis was performed as detailed for annotated transcripts, except differential gene expression was only determined at the transcript level without gene-level summarization. We then searched the differentially expressed transcripts for non-coding RNAs identified in the Rfam database.
Functional classification and gene ontology (GO) enrichment analysis
The Protein Analysis Through Evolutionary Relationships (PANTHER) Classification System  was used to functionally classify the annotated transcriptome. Classification was performed to determine GO slim categories for Biological Process (BP), Molecular Function (MF), Cellular Component (CC). Slim categories were chosen as opposed to traditional GO categories to provide summarized overview of the transcriptome.
Genes that were differentially expressed based both on the gene and transcript-level approaches were analyzed in PANTHER to determine their GO slim BP categories. GO enrichment analysis was also performed on these differentially expressed genes using the Database for Annotation, Visualization and Integrated Discovery (DAVID v6.8) [90, 91]. The annotated transcriptome was used as background. Significant enrichment of GO Biological Process categories was tested and categories were considered enriched if the FDR < 0.05.
Availability of data and materials
The RNA-Seq data generated and analyzed in this study are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA), BioProject ID number PRJNA603387, (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA603387), BioSample accession numbers SAMN13930926-SAMN13930949, Experiment accession numbers SRX7650130-SRX7650153, Run accession numbers SRR10988709-SRR10988686.
Basic local alignment search tool
Database for annotation, visualization and integrated discovery
long non-coding RNAs
Protein analysis through evolutionary relationships
High throughput RNA sequencing
Small nucleolar RNAs, C/D Box
Storey KB, Storey JM. Freeze tolerance and freeze avoidance in Ectotherms. In: Wang LCH, editor. Animal adaptation to cold. Advances in comparative and environmental physiology, vol. 4. Berlin: Springer; 1989. p. 51–82.
Costanzo JP, Lee RE. Avoidance and tolerance of freezing in ectothermic vertebrates. J Exp Biol. 2013;216:1961–7.
Storey KB, Storey JM. Molecular physiology of freeze tolerance in vertebrates. Physiol Rev. 2017;97:623–65.
Storey KB, Storey JM. Physiology, biochemistry and molecular biology of vertebrate freeze tolerance: the wood frog. In: Benson E, Fuller B, Lane N, editors. Life in the frozen state. Boca Raton: CRC Press; 2004. p. 243–74.
Zimmerman SL, Frisbie J, Goldstein DL, West J, Rivera K, Krane CM. Excretion and conservation of glycerol, and expression of aquaporins and glyceroporins, during cold acclimation in Cope's gray tree frog Hyla chrysoscelis. Am J Physiol- Reg I. 2007;292:R544–55.
do MCF A, Lee RE, Costanzo JP. Enzymatic regulation of seasonal glycogen cycling in the freeze-tolerant wood frog, Rana sylvatica. J Comp Physiol B. 2016;186:1045–58.
Mutyam V, Puccetti MV, Frisbie J, Goldstein DL, Krane CM. Dynamic regulation of aquaglyceroporin expression in erythrocyte cultures from cold- and warm-acclimated Cope's gray treefrog, Hyla chrysoscelis. J Exp Zool. 2011;315A:424–37.
Joanisse DR, Storey KB. Oxidative damage and antioxidants in Rana sylvatica, the freeze-tolerant wood frog. Am J Physiol Reg I. 1996;271:R545–53.
Kültz D. Molecular and evolutionary basis of the cellular stress response. Annu Rev Physiol. 2005;67:225–57.
Kiss AJ, Muir TJ, Lee RE Jr, Costanzo JP. Seasonal variation in the hepatoproteome of the dehydration- and freeze-tolerant wood frog, Rana sylvatica. IJMS. 2011;12:8406–14.
Duellman WE, Marion AB, Hedges SB. Phylogenetics, classification, and biogeography of the treefrogs (Amphibia: Anura: Arboranae). Zootaxa. 2016;4104:001–109.
Burkholder G. Hyla chrysoscelis (gray treefrog) hibernacula. Herpetological Rev. 1998;29:231.
Johnson JR, Mahan RD, Semlitsch RD. Seasonal terrestrial microhabitat use by gray treefrogs (Hyla versicolor) in Missouri oak-hickory forests. Herpetologica. 2008;64:259–69.
Costanzo JP, Wright MF, Lee RE. Freeze tolerance as an overwintering adaptation in Cope's grey treefrog (Hyla chrysoscelis). Copeia. 1992;2:565–9.
do MCF A, Frisbie J, Goldstein DL, Krane M. The cryoprotectant system of Cope’s gray treefrog, Dryophytes chrysoscelis: responses to cold acclimation, freezing, and thawing. J Comp Physiol B. 2018;188:611–21.
Irwin JT, Lee RE. Geographic variation in energy storage and physiological responses to freezing in the gray treefrogs Hyla versicolor and H. chrysoscelis. J Exp Biol. 2003;206:2859–67.
Storey JM, Storey KB. Adaptations of metabolism for freeze tolerance in the gray tree frog, Hyla versicolor. Can J Zool. 1985;63:49–54.
Layne JR Jr, Stapleton MG. Annual variation in glycerol mobilization and effect of freeze rigor on post-thaw locomotion in the freeze-tolerant frog Hyla versicolor. J Comp Physiol B. 2008;179:215–21.
Geiss L, MCF DA, Frisbie J, Goldstein DL, Krane CM. Postfreeze viability of erythrocytes from Dryophytes chrysoscelis. J Exp Zool. 2019;331:308–13.
Goldstein DL, Frisbie J, Diller A, Pandey RN, Krane CM. Glycerol uptake by erythrocytes from warm- and cold-acclimated Cope’s gray treefrogs. J Comp Physiol B. 2010;180:1257–65.
Layne JR. Freeze tolerance and cryoprotectant mobilization in the gray treefrog (Hyla versicolor). J Exp Zool. 1999;283:221–5.
Williams JB, Lee RE. Differences in cold tolerance, desiccation resistance, and cryoprotectant production between three populations of Eurosta solidaginis collected from different latitudes. J Comp Physiol B. 2007;178:365–75.
Driedzic WR. Rainbow smelt: the unusual case of cryoprotection by sustained glycerol production in an aquatic animal. J Comp Physiol B. 2015;185:487–99.
Zhang J, Storey KB. Akt signaling and freezing survival in the wood frog, Rana sylvatica. Biochim Biophys Acta. 1830;2013:4828–37.
Keenan SW, Hill CA, Kandoth C, Buck LT, Warren DE. Transcriptomic responses of the heart and brain to anoxia in the western painted turtle. PLoS One. 2015;10:e0131669.
Reilly BD, Schlipalius DI, Cramp RL, Ebert PR, Franklin CE. Frogs and estivation: transcriptional insights into metabolism and cell survival in a natural model of extended muscle disuse. Physiol Genomics. 2013;45:377–88.
Shao Y, Wang L-J, Zhong L, Hong M-L, Chen H-M, Murphy RW, Wu D-D, Zhang Y-P, Che J. Transcriptomes reveal the genetic mechanisms underlying ionic regulatory adaptations to salt in the crab-eating frog. Sci Rep. 2015;5:17551.
The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2017;45:D158–69.
Bryant DM, Johnson K, DiTommaso T, Tickle T, Couger MB, Payzin-Dogru D, Lee TJ, Leigh ND, Kuo TH, Davis FG, Bateman J, Bryant S, Guzikowski AR, Tsai SL, Coyne S, Ye WW, Freeman RM Jr, Peshkin L, Tabin CJ, Regev A, Haas BJ, Whited JL. A tissue-mapped axolotl De novo Transcriptome enables identification of limb regeneration factors. Cell Rep. 2017;18(3):762–76. https://doi.org/10.1016/j.celrep.2016.12.063.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Nawrocki EP, Burge SW, Bateman A, Daub J, Eberhardt RY, Eddy SR, Floden EW, Gardner PP, Jones TA, Tate J, Finn RD. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 2015;43:D130–7.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:31.
Qiao L, Yang W, Fu J, Song Z. Transcriptome profile of the green odorous frog (Odorrana margaretae). PLoS One. 2013;8:e75211.
Zhang W, Guo Y, Li J, Huang L, Kazitsa EG, Wu H. Transcriptome analysis reveals the genetic basis underlying the seasonal development of keratinized nuptial spines in Leptobrachium boringii. BMC Genomics. 2016;17:978.
Huang L, Li J, Anbourkaria H, Luo Z, Zhao M, Wu H. Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Reports. 2016;6:24069. https://doi.org/10.1038/srep24069.
Bansal S, Luu BE, Storey KB. MicroRNA regulation in heart and skeletal muscle over the freeze–thaw cycle in the freeze tolerant wood frog. J Comp Physiol B. 2015;186:229–41.
Hadj-Moussa H, Storey KB. Micromanaging freeze tolerance: the biogenesis and regulation of neuroprotective microRNAs in frozen brains. Cell Mol Life Sci. 2018;75:3635–47.
Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2016;4:1521. https://doi.org/10.12688/f1000research.7563.2.
Zhang D, Tomisato W, Su L, et al. Skin-specific regulation of SREBP processing and lipid biosynthesis by glycerol kinase 5. Proc Natl Acad Sci U S A. 2017;114(26):E5197–206.
Stogsdill B, Frisbie J, Krane CM, Goldstein DL. Expression of the aquaglyceroporin HC-9 in a freeze-tolerant amphibian that accumulates glycerol seasonally. Physiol Rep. 2017;5(15):e13331.
Rosendale AJ, Lee RE Jr, Costanzo JP. Effect of physiological stress on expression of glucose transporter 2 in liver of the wood frog, Rana sylvatica. J Exp Zool. 2014;321:566–76.
Rosendale AJ, Romick-Rosendale LE, Watanabe M, Dunlevy ME, Benoit JB. Mechanistic underpinnings of dehydration stress in the American dog tick revealed through RNA-Seq and metabolomics. J Exp Biol. 2016;219:1808–19.
Storey K. Storey. Heat shock proteins and hypometabolism: adaptive strategy for proteome preservation. RRB. 2011;2:57–68.
Uchida T, Rossignol F, Matthay MA, Mounier R, Couette S, Clottes E, Clerici C. Prolonged hypoxia differentially regulates hypoxia-inducible factor (HIF)-1 and HIF-2 expression in lung epithelial cells: implication of natural antisense HIF-1. J Biol Chem. 2004;279:14871–8.
Kim J, Cha Y-N, Surh Y-J. A protective role of nuclear factor-erythroid 2-related factor-2 (Nrf2) in inflammatory disorders. Mutat Res - Fund Mol M. 2010;690:12–23.
Gracey AY, Fraser EJ, Li W, Fang Y, Taylor RR, Rogers J, Brass A, Cossins AR. Coping with cold: an integrative, multitissue analysis of the transcriptome of a poikilothermic vertebrate. P Natl Acad Sci USA. 2004;101:16970–1691.
Hermes-Lima M, Moreira DC, Rivera-Ingraham GA, Giraud-Billoud M, Genaro-Mattos TC, Campos ÉG. Preparation for oxidative stress under hypoxia and metabolic depression: revisiting the proposal two decades later. Free Radical Bio Med. 2015;89:1122–43.
Elmore S. Apoptosis: a review of programmed cell death. Toxicol Pathol. 2007;35:495–516.
Owens DM, Keyse SM. Differential regulation of MAP kinase signaling by dual-specificity protein phosphatases. Oncogene. 2007;26:3203–13.
Salvador JM, Brown-Clay JD, Fornace AJ. Gadd45 in stress signaling, cell cycle control, and apoptosis. In: Advances in experimental medicine and biology. New York: Springer; 2013. p. 1–19.
Cheng C, Alexander R, Min R, Leng J, Yip KY, Rozowsky J, Yan KK, Dong X, Djebali S, Ruan Y, Davis CA, Carninci P, Lassman T, Gingeras TR, Guigo R, Birney E, Weng Z, Snyder M, Gerstein M. Understanding transcriptional regulation by integrative analysis of transcription factor binding data. Genome Res. 2012;22:1658–67.
Bluyssen HA, Levy DE. Stat2 is a transcriptional activator that requires sequence-specific contacts provided by stat1 and p48 for stable interaction with DNA. J Biol Chem. 1997;272:4600–5.
Nagasawa K, Tanizaki Y, Okui T, Watarai A, Ueda S, Kato T. Significant modulation of the hepatic proteome induced by exposure to low temperature in Xenopus laevis. Biol Open. 2013;2:1057–69.
Richards RC, Short CE, Driedzic WR, Ewart KV. Seasonal changes in hepatic gene expression reveal modulation of multiple processes in rainbow smelt (Osmerus mordax). Mar Biotechnol. 2010;12:650–63.
Todgham AE, Crombie TA, Hofmann GE. The effect of temperature adaptation on the ubiquitin–proteasome pathway in notothenioid fishes. J Exp Biol. 2017;220:369–78.
Lynch M, Marinov GK. The bioenergetic costs of a gene. Proc Natl Acad Sci U S A. 2015;112:15690–5.
Yu F, Deng H, Yao H, Liu Q, Su F, Song E. Mir-30 reduction maintains self-renewal and inhibits apoptosis in breast tumor-initiating cells. Oncogene. 2010;29:4194–204.
Tu X, Zheng X, Li H, Cao Z, Chang H, Luan S, Zhu J, Chen J, Zang Y, Zhang J. MicroRNA-30 protects against carbon tetrachloride-induced liver fibrosis by attenuating transforming growth factor beta signaling in hepatic stellate cells. Toxicol Sci. 2015;146:157–69.
Scott MS, Ono M. From snoRNA to miRNA: dual function regulatory non-coding RNAs. Biochimie. 2011;93:1987–92.
Merchen TD, Boesen EI, Gardner JR, Harbarger R, Kitamura E, Mellor A, Pollock DM, Ghaffari A, Podolsky R, Nahman NS Jr. Indoleamine 2,3-dioxygenase inhibition alters the non-coding RNA transcriptome following renal ischemia-reperfusion injury. Transpl Immunol. 2014;30:140–4.
Quinn JJ, Chang HY. Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet. 2016;217(1):47–62. https://doi.org/10.1038/nrg.2015.10 Review.
Wilusz JE. Long noncoding RNAs: re-writing dogmas of RNA processing and stability. Biochim Biophys Acta. 2016;1859(1):128–38. https://doi.org/10.1016/j.bbagrm.2015.06.003.
Hammond SA, Warren RL, Vandervalk BP, Kucuk E, Khan H, Gibb EA, Pandoh P, Kirk H, Zhao Y, Jones M, Mungall AJ, Coope R, Pleasance S, Moore RA, Holt RA, Round JM, Ohora S, Walle BV, Veldhoen N, Helbing CC, Birol I. The north American bullfrog draft genome provides insight into hormonal regulation of long noncoding RNA. Nat Commun. 2017;8(1):1433. https://doi.org/10.1038/s41467-017-01316-7.
Calixto CPG, Tzioutziou NA, James AB, Hornyik C, Guo W, Zhang R, Nimmo HG, Brown JWS. Cold-dependent expression and alternative splicing of Arabidopsis long non-coding RNAs. Front Plant Sci. 2018;10:235. https://doi.org/10.3389/fpls.2019.00235.
Kindgren P, Ard R, Ivanov M, Marquardt S. Transcriptional read-through of the long non-coding RNA SVALKA governs plant cold acclimation. Nat Commun. 2018;9(1):4561. https://doi.org/10.1038/s41467-018-07010-6.
Ran MX, Li Y, Zhang Y, Liang K, Ren YN, Zhang M, Zhou GB, Zhou YM, Wu K, Wang CD, Huang Y, Luo B, Qazi IH, Zhang HM, Zeng CJ. Transcriptome sequencing reveals the differentially expressed lncRNAs and mRNAs involved in cryoinjuries in frozen-thawed Giant panda (Ailuropoda melanoleuca) sperm. Int J Mol Sci. 2018;19(10):3066. https://doi.org/10.3390/ijms19103066.
Jansen HT, Trojahn S, Saxton MW, et al. Hibernation induces widespread transcriptional remodeling in metabolic tissues of the grizzly bear. Commun Biol. 2019;2:336. https://doi.org/10.1038/s42003-019-0574-4.
Vermillion KL, Anderson KJ, Hampton M, Andrews M. Gene expression changes controlling distinct adaptations in the heart and skeletal muscle of a hibernating mammal. Physiol Genomics. 2015;47(3):58–74.
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(8):2413–26. https://doi.org/10.1093/gbe/evw163.
Metzger DCH, Schulte PM. Persistent and plastic effects of temperature on DNA methylation across the genome of threespine stickleback (Gasterosteus aculeatus). Proc Biol Sci. 2017;284(1864):20171667. https://doi.org/10.1098/rspb.2017.1667.
Nitzan T, Kokou F, Doron-Faigenboim A, Slosman T, Biran J, Mizrahi I, Zak T, Benet A, Cnaani A. Transcriptome analysis reveals common and differential response to low temperature exposure between tolerant and sensitive blue Tilapia (Oreochromis aureus). Front Genet. 2019;10:100. https://doi.org/10.3389/fgene.2019.00100.
MacMillan HA, Knee JM, Dennis AB, Udaka H, Marshall KE, Merritt TJ, Sinclair BJ. Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome. Sci Rep. 2016;6:28999. https://doi.org/10.1038/srep28999.
Lou X, Wang H, Ni X, Gao Z, Iqbal S. Integrating proteomic and transcriptomic analyses of loquat (Eriobotrya japonica Lindl.) in response to cold stress. Gene. 2018;677:57–65. https://doi.org/10.1016/j.gene.2018.07.022.
Wang K, Bai ZY, Liang QY, Liu QL, Zhang L, Pan YZ, Liu GL, Jiang BB, Zhang F, Jia Y. Transcriptome analysis of chrysanthemum (Dendranthema grandiflorum) in response to low temperature stress. BMC Genomics. 2018;19(1):319. https://doi.org/10.1186/s12864-018-4706-x.
Hannah MA, Heyer AG, Hincha DK. A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genet. 2005;1(2):e26.
de Boer TE, Roelofs D, Vooijs R, Holmstrup M, Amorim MJB. Population-specific transcriptional differences associated with freeze tolerance in a terrestrial worm. Ecol Evol. 2018;8(7):3774–86. https://doi.org/10.1002/ece3.3602.
Toxopeus J, Des Marteaux LE, Sinclair BJ. How crickets become freeze tolerant: the transcriptomic underpinnings of acclimation in Gryllus veletis. Comp Biochem Physiol Part D Genomics Proteomics. 2019;29:55–66. https://doi.org/10.1016/j.cbd.2018.10.007.
Parker DJ, Vesala L, Ritchie MG, Laiho A, Hoikkala A, Kankare M. How consistent are the transcriptome changes associated with cold acclimation in two species of the Drosophila virilis group? Heredity. 2015;115(1):13–21. https://doi.org/10.1038/hdy.2015.6.
Zuther E, Schaarschmidt S, Fischer A, Erban A, Pagter M, Mubeen U, Giavalisco P, Kopka J, Sprenger H, Hincha DK. Molecular signatures associated with increased freezing tolerance due to low temperature memory in Arabidopsis. Plant Cell Environ. 2019;42(3):854–73. https://doi.org/10.1111/pce.13502.
Aguilar OA, Hadj-Moussa H, Storey KB. Freeze-response regulation of MEF2 proteins and downstream gene networks in muscles of the wood frog, Rana sylvatica. J Thermal Biol. 2017;67:1–8.
Wu CW, Tessier SN, Storey KB. Stress-induced antioxidant defense and protein chaperone response in the freeze-tolerant wood frog Rana sylvatica. Cell Stress Chaperones. 2018;23:1205–17.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, MacManes MD, Ott M, Orvis J, Pochet N, Strozzi F, Weeks N, Westerman R, William T, Dewey CN, Henschel R, LeDuc RD, Friedman N, Regev A. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Bray N, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7.
Team RC. R: A language and environment for statistical computing. [Online]. https://www.R-project.org/.
Cook RD. Detection of influential observation in linear regression. Technometrics. 1977;19:15–8.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc B Met. 1995;57:289–300.
Mi H, Poudel S, Muruganujan A, Casagrande JT, Thomas PD. PANTHER version 10: expanded protein families and functions, and analysis tools. Nucleic Acids Res. 2016;44:D336–42.
Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.
We thank Dr. Andrew Rosendale for advice regarding transcriptome annotation and for providing constructive comments on an earlier version of the manuscript. Current affiliation for R.C.: University of Pittsburgh School of Medicine, Medical Scientist Training Program, 3550 Terrace St, Pittsburgh, PA 15261.
This work was supported by NSF IOS-1121457 to D.L.G. and C.K., and Schuellein Chair in the Biological Sciences to C.K., and Schuellein Chair Postdoctoral Fellowship to M.C.F.A. These funding bodies had no role in the design of the study, collection, analysis, and interpretation of data, or in writing the manuscript.
Ethics approval and consent to participate
Rearing and all experiments were conducted with the approval of the Laboratory Animal Care and Use Committee at Wright State University.
Consent for publication
The authors declare they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
genes differentially regulated in Dryophytes chrysoscelis compared to warm animals with corresponding gene symbol, UniprotID, Transcript ID, and fold-change in each treatment. Table S2. genes differentially expressed in frozen Dryophytes chrysoscelis, relative to cold animals, with corresponding gene symbol, UniprotID, Transcript ID, and fold-change.
About this article
Cite this article
do Amaral, M.C.F., Frisbie, J., Crum, R.J. et al. Hepatic transcriptome of the freeze-tolerant Cope’s gray treefrog, Dryophytes chrysoscelis: responses to cold acclimation and freezing. BMC Genomics 21, 226 (2020). https://doi.org/10.1186/s12864-020-6602-4