1 Citation   Views   Downloads

Navigating oxygen deprivation: liver transcriptomic responses of the red eared slider turtle to environmental anoxia

View article
RT @jzhangcad: My article “Navigating oxygen deprivation: liver transcriptomic responses of the red eared slider turtle to environmental an…
RT @jzhangcad: My article “Navigating oxygen deprivation: liver transcriptomic responses of the red eared slider turtle to environmental an…
Navigating oxygen deprivation: liver transcriptomic responses of the red eared slider turtle to environmental anoxia https://t.co/QKpUq8DpFu @thePeerJ https://t.co/dCP6TWs0Gd
90 days ago
My article “Navigating oxygen deprivation: liver transcriptomic responses of the red eared slider turtle to environmental anoxia” has been published today in @thePeerJ https://t.co/jUKWgq6jTU #Biochemistry #MolecularBiology #Zoology
Biochemistry, Biophysics and Molecular Biology


Maintaining optimal and continuous oxygen supply to organs (in order to maintain high rates of aerobic metabolism) is not universal across the animal kingdom. Many facultative anaerobic organisms, including a variety of ectothermic vertebrates, are able to survive extended periods without a supply of oxygen (Biggar, Groom & Storey, 2012; Storey, 2007; Nilsson & Renshaw, 2004). For example, a lack of available oxygen is common among turtles, which spend the majority of their lives diving underwater (Clark & Miller, 1973; Jackson, 1968). For bouts of short-term oxygen deprivation, such as during underwater dives for food, simply increasing anaerobic metabolism is able to meet metabolic demands (Storey & Storey, 2004). However, for long term survival, such as overwintering in ice-locked ponds and waterways, turtles must resort to at least one of two established modes of survival: extrapulmonary modes of oxygen uptake, or a severe reduction in metabolic demand such that anaerobic metabolism can meet energetic demand (i.e., hypometabolism).

The best facultative anaerobes among vertebrates are members of the genera Trachemys (pond slider turtles) and Chrysemys (painted turtles). These turtles are able to survive without oxygen for up to two weeks at ∼16 °C and for 12 to 18 weeks at ∼3 °C (Ultsch, 1985; Ultsch & Jackson, 1982). By comparison, mammalian tissues are notoriously sensitive to even brief episodes of anoxia; for example, humans can endure only a few minutes of anoxia without irreversible brain damage. Previously reported by Hochachka et al. (1996), the main features of anoxia survival via hypometabolism include: (1) a mechanism of oxygen sensing such that anoxic transitions can be signaled to cells, (2) a set of pro-survival genes that are up-regulated to help coordinate adaptation to the hypometabolic state, (3) a set of non-essential genes and (4) non-essential pathways that are down-regulated to conserve energy stores, (5) reduction in activity and impulse frequency in neural tissues, and (6) a sustainable re-balance of metabolism such that cellular energy demands can be met by anaerobic energy production (Hochachka, Land & Buck, 1997). Collectively, these are thought to be fundamental, and highly regulated, processes that allow turtles to survive extended periods of reduced oxygen in a hypometabolic state. In turtles, a profound metabolic rate depression to only 10–20% of the corresponding aerobic resting rate, at the same temperature, occurs in response to anoxia (Storey & Storey, 1990). Indeed, to help establish this metabolically depressed state, it has been well established that turtles respond to anoxia with a suit of cellular modifications that involve transcriptional, translational, post-translational (i.e., reversible protein phosphorylation), as well as epigenetic, and microRNA-dependent post-transcriptional regulations (Storey, 2007; Krivoruchko & Storey, 2015; Biggar & Storey, 2017).

Upon encountering hypoxic conditions, the first response of the turtle is to enhance oxygen delivery and the oxygen extraction systems. If oxygen concentrations continue to fall, the systemic alterations to oxygen extraction quickly become inadequate to supply enough oxygen to deprived tissues. When this occurs, oxygen-independent metabolic pathways, such as anaerobic glycolysis, are fully recruited and are followed by regulatory actions that reduce oxygen demand by cells (Hochachka et al., 1996). This introduces a very important issue—increasing the rate of glycolysis does increase ATP output, but also results in a quick depletion of internal fuel reserves, as well as a large accumulation of acidic end products that must be dealt with (Kelly & Storey, 1988). One predominant component of ATP conservation to match the supply provided by anaerobic glycosis is the suppression of protein synthesis.

Protein synthesis consumes a very large portion of available ATP turnover in cells under normoxic conditions, requiring about 5 ATP equivalents per peptide bond of a newly synthesized protein (Cramer et al., 1991). Some freshwater turtles have been shown to decrease the rate of ATP use by suppressing protein synthesis to as low as 6% during anoxia (Hochachka et al., 1996). Several studies have explored the in vivo protein synthesis rates during anoxia-induced metabolic depression in turtles (Brooks & Storey, 1993; Land & Hochachka, 1994; Fraser et al., 2001). Such a dramatic regulation of the transcriptome during hypometabolism might be achieved in several ways: (1) a reduction of the amount of mRNA transcripts available (i.e., changes in mRNA expression), (2) differential regulation of the assembly and activity of ribosomal translational machinery, (3) epigenetic-level gene regulation, and (4) microRNA interference.

To date, a large body of literature has established that in response to anoxia, the red-eared slider turtle is able to respond to oxygen lack through a number of molecular mechanisms that include; (1) the expression of various molecular chaperones, (2) increased antioxidant defense, (3) stress-responsive transcription factors, (4) pro-survival proteins, and (5) a significant restructuring of metabolism such that ATP demand can be met in the new anaerobic environment. For example, one recent study documented epigenetic changes that dynamically occur on histone proteins to alter the chromatin state in a manner that is responsive to anoxia. This suggests chromatin-level regulation of gene expression in an anoxia-responsive manner in the turtle (Wijenayake, Hawkins & Storey, 2018). Given the large amount of supporting knowledge behind the molecular mechanisms that permits this turtle to survive extended periods of anoxia, we chose the red-eared slider turtle as the most appropriate model for this study. Supplemented with the recent sequencing of the western painted turtle genome (C. picta bellii), an unprecedented degree of insight can now be gained into the transcriptomic response behind turtle anoxia tolerance in general (Shaffer et al., 2013; Keenan et al., 2015; Hong et al., 2019). In the present study, we utilize this genomic resource for several purposes that include: (1) to investigate whether anoxic-responsive changes in the transcriptome are reflective of known physiological responses, (2) to provide insight into transcriptome-wide changes in biological processes and molecular pathways that are sensitive to anoxia, and (3) to identify new targets for future functional study.

Material and Methods

Animal care and treatment

Adult female red-eared sliders (T. s. elegans), 700–1,500 g, were acquired from local suppliers (Ottawa, Ontario, Canada) and held at 5 ± 1 °C in large 50 L plastic tanks (2 turtles per tank) filled with dechlorinated tap water for an average of 7–8 days before use. Control (normoxic) turtles were sampled from this condition. For anoxia exposure, turtles were transferred to large buckets at 5 ± 1 °C that had been previously bubbled with N2 gas for 1 hr; 2 turtles were added per bucket in 30 min intervals to allow for sufficient sampling and post-processing time between individuals. Bubbling was continued for 1 hr after the last turtle was added and was reinitiated again during sampling of the animals. A wire mesh was fitted into the tank about 5 cm below the water surface so that turtles remained submerged throughout the 20 hrs experimental anoxia exposure. All animals were killed by decapitation and liver tissue was rapidly dissected out, frozen in liquid nitrogen and stored at −80 °C until use. The time-course of anoxia exposure was chosen as described in Jackson (2002). All animals were cared for in accordance with the guidelines of the Canadian Council on Animal Care and all experimental procedures had the prior approval of the Carleton University Animal Care Committee (IACUC approval #B09-20).

RNA isolation and library preparation

Total RNA was isolated from the liver tissue of control and 20 hrs anoxic turtles (n = 3 individual turtles) using Trizol (Invitrogen; Cat#15596-018). Briefly, 100 mg of tissue was homogenized in 1 mL Trizol using a Polytron homogenizer followed by the addition of 200 µL of chloroform and centrifugation at 10,000×g for 15 min at 4 °C. The upper aqueous layer (containing RNA) was removed and collected into an RNAse-free microcentrifuge tube. Total RNA was then precipitated with the addition of 500 µL of isopropanol followed by incubation for 10 min at room temperature. Samples were then centrifuged again as above. The RNA pellet was washed with 70% ethanol. Samples were loaded onto Qiagen RNeasy filter columns, and kit protocol was followed according to manufacturer’s instructions, including an on-column DNase step to remove any residual contaminating DNA. RNA integrity and concentration were assessed using the Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) at the McGill University and Genome Quebec Innovation Centre (Montreal, QC). RNA Integrity Number (RIN) values for total RNA used for sequencing were greater than 7. Libraries were constructed with the TruSeq Stranded mRNA Sample Preparation Kit (Illumina) following manufacturers guidelines.

Transcriptome sequencing and assembly

The RNA sequencing performed at the McGill University and Genome Quebec Innovation Centre (Montreal, QC) generated 148 to 178 million paired reads per library using a Illumina HiSeq 2500 T1 sequencer. Trimming and clipping of Illumina sequencing adapters and other adapters was performed by the sequencing facility using Trimmomatic software (Bolger, Lohse & Usadel, 2014). Reads are trimmed from the 3′end to have a phred score of at least 30. Illumina sequencing adapters are removed from the reads, and all reads are required to have a length of at least 32. The filtered reads were aligned to the Chrysemys_picta_bellii-3.0.3 reference genome (NCBI) using TopHat software (Trapnell, Pachter & Salzberg, 2009). TopHat splits reads to align them across known and novel splice junctions. The filtered and aligned T. scripta elegans RNAseq dataset represents a minimum 49x coverage of the C. picta bellii exome. To estimate transcript and gene abundances, Cufflinks software (v1.3.0) was used (Trapnell et al., 2013), generating normalized FPKMs (Fragments Per Kilobase of exon model per Million mapped fragments) for each annotated gene. All trimming and alignment statistics are provided in Table 1. All raw sequencing reads are able available through the NCBI sequence read archive database (SRA PRJNA525986). The GFF used for alignment is available at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA78657/.

Table 1:
Trimming and alignment statistics.
Alternative alignments indicates the number of duplicate read entries providing alternative coordinates.
Sample Raw reads Surviving reads % Aligned reads % Alternative alignments % Coverage Exonic rate Genes
C-1 157,384,476 153,070,220 97 96,895,249 63 11,227,181 12 63 0.73 19,230
C-2 177,029,988 171,370,576 97 101,933,083 59 10,532,358 10 63 0.75 18,500
C-3 160,902,862 155,600,064 97 91,604,554 59 9,145,947 10 57 0.75 18,309
20-1 161,001,556 155,718,808 97 85,453,907 55 8,975,060 11 54 0.74 18,178
20-2 148,150,044 143,300,486 97 78,889,776 55 8,167,035 10 51 0.74 18,169
20-3 150,765,402 145,901,344 97 80,739,759 55 8,289,079 10 49 0.74 18,252
DOI: 10.7717/peerj.8144/table-1

Gene differential expression (DE) analysis

Normalized expression values were subject to unsupervised hierarchical cluster analysis. A heat map was used to visualize the clustering results and the overall gene expression pattern for each experimental condition and individual animal sample. All raw FPKM values were transformed through log10(FPKM+1) prior to plotting. The R package gplots (Warnes et al., 2016) with a custom workflow was used for generating the heat map. Volcano distribution was used to visualize the DE (i.e., individual gene level) results. The log2 transformed fold changes were plotted against –log10 transformed p-values (false discovery rate, or FDR, adjusted) for each gene.

Gene set (GS) analysis

Gene Ontology (GO) term and KEGG pathway (Ashburner et al., 2000; Kanehisa et al., 2016) were used for GS analysis. In order to explore the molecular events responsible for the potential anoxia-responsive adjustments to the cell, both biological processes (BP) and molecular function (MF) sets were used for GO term analysis. KEGG gene sets were used to further explore the pathways affected by the anoxic conditions and generate visual representations with the integration of the relevant DE results.

Gene level statistics of all genes were used for GS analysis, as opposed to only significant genes (i.e., the contingency table method). As such, the current method was able to take the overall gene expression profile into consideration for GS enrichment without relying on the univariate analysis significance threshold. By combining gene level p-values with directional statistics (e.g., fold change or t-statistics), the present GS assessment provided both directional and non-directional GS statistics. Specifically, non-directional, distinct directional and mixed directional tests were conducted on both GO term and KEGG GS analyses. With the distinct-directional analysis, we were able to explore the overall anoxia response of the pathways, as the test weighs the opposite directional changes exhibited by subsets of the genes, thereby determining the direction in which the gene sets are regulated. The mixed directional test investigates the directional impact to the pathways from the subsets of the genes included in the GS based on DE statistics. For example, a pathway can be regulated in both positive and negative directions if the two opposite directional regulations on two subsets of genes significantly affect the GS statistics. For both the distinct and mixed directional tests, GS level p-values were calculated based on the directionality exhibited in the gene level statistics. The non-directional GS analysis identified gene sets that were significantly impacted by the anoxia stress solely based on the significance test results (p-values) from the DE results.

To minimize analysis bias, several GS enrichment methods were used to comprehensively and thoroughly assess the biological significance of anoxia-responsive gene expression profiles through a consensus score-based ranking method (Väremo, Nielsen & Nookaew, 2013). A total of eight enrichment methods were used, including parametric (Fisher’s combined probability, Stouffer’s test, Reporter features test, and parametric analysis of gene set enrichment (i.e., PAGE) test) (Patil & Nielsen, 2005; Kim & Volsky, 2005), non-parametric (Tail strength test, Wilcoxon rank-sum test, and gene set enrichment analysis (i.e., GSEA) test) approaches (Taylor & Tibshirani, 2006; Smyth, 2005; Subramanian et al., 2005), as well as a maxmean method developed by Efron and Tibshirani (Efron & Tibshirani, 2007). Gene sampling permutation test (x1000) was used for generating the null distribution for all enrichment methods to determine the significance of the enrichment results (Goeman & Buhlmann, 2007). All the GS p-values from the individual enrichment methods were adjusted by FDR. Consistent with the DE analysis, best p-value was set for the gene level statistics for repeated genes. All the gene sets were ranked for individual enrichment tests using the associated p-values. The median value of the ranks was considered the consensus score (normalized to a 0–100 scale) for each GS.

Anoxia-induced transcriptomic changes.

Figure 1: Anoxia-induced transcriptomic changes.

(A) Heat map showing both the results for unsupervised hierarchical clustering analysis and overall gene expression pattern. FPKM values were plotted after log10(FPKM+1) transformation. (B) A volcano plot showing the results for DE analysis. Logarithmic transformed p-values (FDR corrected) and fold changes are plotted. Orange colour represents the fold changes equal or greater than 1.5; red represents p-values less than 0.05; green represents the genes that meet both criteria.

Quadrant scatterplot and boxplot were used to visualize the results for GS analysis for all five directionality classes. Specifically, the median GS p-value (cut-off: 0.1) were calculated based on the multiple GS enrichment tests and plotted against the GS consensus score in the boxplots, and used as the consensus and final p-value for the gene sets tested. Gene sets with a final p-value less than 0.05 were considered significantly altered by anoxia. Additionally, gene sets from either GO term or KEGG pathways were sorted according to the consensus score (cut-off: 50) and plotted using boxplot. Furthermore, selected pathways were visualized with the related DE results. The R package pathview (Luo & Brouwer, 2013) was used for ID conversion and KEGG pathway visualization. GS enrichment tests and consensus score processing was conducted using the R package piano (Väremo, Nielsen & Nookaew, 2013) with a custom workflow specific to the current study. Results visualization for GS analysis was carried out using R packages reshape2 (Wickham, 2007) and ggplot2 (Wickham, 2009) with customization.


Gene differential expression (DE) analysis

To investigate the effects of anoxia on the transcriptomic responses of the liver, adult turtles were submerged in anoxic water for 20 hrs at 5 °C, and the liver tissue was sampled and analyzed using RNA-seq. A heat map was used to visualize unsupervised hierarchical clustering and overall gene expression profiles (Fig. 1A). While variance existed between individual samples in the same experimental group, the clustering analysis demonstrated a clear separation of the overall gene expression pattern between control and anoxic turtles. A volcano plot was used to visualize DE analysis results (Fig. 1B). A total of 3,105 genes showed statistically significant changes in gene expression between the two animal groups, of which 971 also exhibited a fold change equal to or greater than 50% of control normoxic values. The plot uses red colour for genes with an adjusted p-value less than 0.05 (statistically significant), orange for fold changes greater than 50% of control normoxic values, and green for the genes that meet both thresholds. Within this dataset, expression of 64 genes increased greater than 5-fold when compared to control values (Table 2), while 417 genes decreased to less than 20% of control values. All fold changes in gene expression are listed as Table S1.

Table 2:
FPKM values measured for genes found to increase greater than 5-fold in response to anoxia.
Fold change Gene name Gene ID Locus Mean FPKM (Control) Mean FPKM (Anoxia)
5.01 Heterogeneous nuclear ribonucleoprotein C gene16690 NW_007281930.1:7005-18045 1.00 4.64
5.12 Cholesteryl ester transfer protein gene3179 NC_024229.1:1268894-1284850 1.51 7.48
5.25 Condensin-2 complex subunit H2 gene12017 NW_007281530.1:1767412-1786677 1.53 7.88
5.25 V-type proton ATPase subunit S1-like protein gene3838 NW_007281339.1:7340969-7369713 0.21 1.03
5.28 Centrosomal protein of 131 kDa gene4049 NW_007281340.1:11239101-11269228 0.43 2.07
5.28 Phosphoserine phosphatase gene20113 NW_007359849.1:8379950-8401297 0.86 4.40
5.30 Tubulin polyglutamylase TTLL6 gene21468 NW_007359872.1:2397248-2424904 0.05 0.25
5.30 RNA helicase Mov10l1 gene1061 NC_024219.1:57385859-57457741 0.11 0.64
5.34 DNA (cytosine-5)-methyltransferase 3-like gene7912 NW_007281399.1:3261284-3279562 0.06 0.34
5.41 Mitochondrial fission regulator 2 gene1488 NC_024220.1:49284333-49299028 0.87 5.47
5.42 Caspase recruitment domain-containing protein 9 gene6122 NW_007281364.1:6381680-6402673 0.23 1.12
5.52 Ectodysplasin-A receptor-associated adapter protein gene6427 NW_007281370.1:1143933-1158884 0.36 1.80
5.53 Apoptosis-inducing, TAF9-like domain 1 gene3583 NC_024234.1:6047045-6056630 0.77 4.12
5.63 NADPH oxidase 1 gene2721 NC_024225.1:13666090-13681627 0.07 0.39
5.70 Occludin gene3858 NW_007281339.1:8859408-8890164 0.85 4.48
5.73 Huntingtin-interacting protein 1 gene12854 NW_007281570.1:421065-557723 0.55 3.36
5.75 Mitochondrial folate transporter/carrier gene10067 NW_007281454.1:2041943-2060033 0.56 3.19
5.83 Exportin-T gene8389 NW_007281409.1:517438-594831 1.21 6.65
5.91 Probable C-mannosyltransferase DPY19L1 gene10618 NW_007281477.1:1919781-2027394 0.64 3.56
5.98 Ras-related protein Rab-26 gene23555 NW_007359899.1:14273841-14485260 0.82 4.25
6.02 Transmembrane protein 218 gene8138 NW_007281402.1:4074561-4089102 0.38 2.05
6.06 Cilia- and flagella-associated protein 251 gene22360 NW_007359884.1:8972215-9006606 0.06 0.39
6.15 Transient receptor potential cation channel subfamily V member 1 gene5628 NW_007281358.1:4608333-4641788 0.14 0.87
6.22 T-box transcription factor TBX6 gene15903 NW_007281789.1:246506-252000 0.33 1.98
6.31 Prominin-1 gene11980 NW_007281529.1:386175-473658 0.10 0.64
6.39 Breast carcinoma-amplified sequence 1 gene9931 NW_007281448.1:1292657-1386845 0.10 0.69
6.40 HERV-H LTR-associating protein 1 gene23918 NW_007359902.1:15300963-15351050 0.14 0.91
6.42 TBC1 domain family member 16 gene4064 NW_007281340.1:12087218-12134952 1.63 9.92
6.47 Sulfotransferase gene22400 NW_007359884.1:10968341-11116484 0.25 1.54
6.81 Putative sodium-coupled neutral amino acid transporter 8 gene3200 NC_024229.1:2867872-2904393 0.07 0.50
6.91 Polypeptide N-acetylgalactosaminyltransferase gene19712 NW_007359841.1:19853527-20196374 0.07 0.48
6.96 Mitochondrial enolase superfamily member 1 gene19866 NW_007359845.1:5442404-5483366 1.96 13.11
7.08 Apical endosomal glycoprotein gene11563 NW_007281514.1:751300-805295 0.04 0.27
7.15 Homogentisate 1,2-dioxygenase gene108 NC_024218.1:16807818-16855000 0.26 1.82
7.24 Ankyrin repeat domain-containing protein 50 gene6263 NW_007281367.1:20195-70713 0.42 2.96
7.29 Hemicentin-2; Receptor protein-tyrosine kinase gene5857 NW_007281360.1:4685366-4806403 0.53 3.90
7.30 Oncostatin-M-specific receptor subunit beta gene5251 NW_007281354.1:2634377-2714975 0.86 6.38
7.39 Intermediate conductance calcium-activated potassium channel protein 4 gene15763 NW_007281776.1:338903-352616 0.16 1.22
7.43 E3 ubiquitin-protein ligase MSL2 gene11846 NW_007281524.1:1835898-1890751 0.29 1.86
7.76 NHS-like protein 2 gene13609 NW_007281610.1:13311-187817 0.47 3.20
8.35 Tetratricopeptide repeat protein 31 gene12423 NW_007281546.1:900243-921746 0.97 7.87
8.35 Transmembrane protein 240 gene3643 NC_024234.1:10334661-10371762 0.12 0.98
9.11 GTP-binding protein Di-Ras3 gene2948 NC_024225.1:40419019-40426191 0.04 0.32
9.21 Ectonucleotide pyrophosphatase/phosphodiesterase family member 1 gene1524 NC_024220.1:52582788-52659285 0.05 0.40
9.45 Integrin alpha-E gene5634 NW_007281358.1:4706392-4771983 0.45 4.04
9.60 Regulator of G-protein signaling 9-binding protein gene3306 NC_024229.1:16667911-16670807 0.03 0.29
9.62 Dpy-19 Like C-Mannosyltransferase 1 gene10618 NW_007281477.1:1919781-2027394 1.46 14.31
9.93 Potassium voltage-gated channel subfamily V member 2 gene2222 NC_024223.1:12230785-12238276 0.05 0.47
10.41 Serine/threonine-protein kinase NIM1 gene5286 NW_007281354.1:7715828-7732336 0.07 0.76
11.11 Fibronectin type III domain-containing protein 4 gene15588 NW_007281760.1:51390-91029 2.19 23.79
11.68 Estradiol 17-beta-dehydrogenase 2 gene5112 NW_007281352.1:6859463-6870938 0.07 0.80
11.75 Pannexin gene8126 NW_007281402.1:2761308-2777293 0.05 0.61
12.25 Nebulin gene8448 NW_007281410.1:3153891-3365730 0.09 1.07
12.42 Perilipin gene5450 NW_007281356.1:5936040-5951851 0.32 3.73
12.82 Leucine-rich repeat-containing protein 38 gene3545 NC_024234.1:4003997-4030353 0.03 0.34
16.22 Tyrosine-protein kinase; Tyrosine-protein kinase ABL2 gene2786 NC_024225.1:22587197-22668610 0.54 8.08
17.04 Glucokinase regulatory protein gene17516 NW_007282306.1:7094-42264 0.64 10.68
18.09 Tripartite motif-containing protein 16; E3 ubiquitin-protein ligase TRIM16 gene22614 NW_007359887.1:6894167-6907028 0.02 0.38
19.71 EF-hand calcium-binding domain-containing protein 3 gene13683 NW_007281615.1:460307-476141 0.04 0.69
22.70 Guanylyl cyclase-activating protein 2 gene8849 NW_007281419.1:2399418-2408781 0.14 3.06
24.76 Contactin-associated protein-like 2 gene20384 NW_007359856.1:12726798-14405054 0.05 1.07
28.83 F-box only protein 2 gene3570 NC_024234.1:4909708-4916629 0.05 1.46
46.44 WD repeat-containing protein 27 gene1362 NC_024220.1:30833693-31040151 0.04 2.05
55.88 Transmembrane protein 198 gene14453 NW_007281670.1:215832-222125 0.05 2.77
DOI: 10.7717/peerj.8144/table-2

Gene-set analysis

To determine the anoxia-responsive regulation of distinct pathways identified from transcriptome data of control and 20 hr anoxic sample groups, we employed a GS analysis. This analysis was used to provide an unbiased insight into common pathways under anoxia-responsive transcriptional control that may aid in low oxygen adaptation in the turtle. Scatterplots were used to display the overall GS results for all directionality classes (non-directional, distinct directional and mixed directional) for both GO terms and KEGG pathways (Fig. 2). The upper right quadrant includes the pathways with a median p-value less than 0.05 and a consensus score equal or less than 50. As shown in Fig. 2A, a total of 246 median p-values (within the five directionality classes) were revealed representing the top significantly enriched GO biological process (BP) terms. Specifically, 49 significantly enriched GO BP terms had a consensus score equal or less than 50 for distinct down-regulation, mixed up-regulation and non-directional regulation classes, while 48 and 51 GO BP terms were enriched for distinct up-regulation and mixed down-regulation classes, respectively (Table S2). Regarding the GO molecular function (MF) terms, 218 median p-values were presented in the upper right quadrant (Fig. 2B), including 52 for mixed down-regulation and non-directional classes, 36 for distinct up-regulation, 47 for distinct down-regulation, and 31 for mixed up-regulation (Table S3). Moreover, KEGG pathway analysis showed 121 median p-values that were less than 0.05 for the top enriched pathways (Fig. 2C); this included 13, 48, 17, 21, and 22 for distinct up-regulation, distinct down-regulation, mixed up-regulation, mixed down-regulation, and non-directional classes, respectively (Table S4). Figure 3 shows boxplots of the top ranked gene sets from both the up and down-directional classes for KEGG pathway enrichment. Boxplots for the rest of the directionality classes can be viewed in the supplemental material (Figs. S1S8).

Quadrant scatterplots showing GS analysis results for five directionality classes.

Figure 2: Quadrant scatterplots showing GS analysis results for five directionality classes.

The upper right quadrant contains the gene sets with a GS median p-value less than 0.05 and a consensus score equal to or less than 50. The total number of the gene sets presented in the upper right quadrant is indicated. (A) GO BP term, (B) GO MF term, and (C) KEGG pathway.
Boxplots showing gene set ranking.

Figure 3: Boxplots showing gene set ranking.

Gene sets are ranked based on the GS consensus score and sorted in ascending order (y-axis, top down). The box representing the individual ranks from GS analyses with the median value (i.e., GS consensus score) displayed as a vertical line. The dots are outliers. Data are for anoxia induced (A) up-regulated distinct directional and (B) down-regulated distinct directional KEGG pathways.

In this study, 25 KEGG pathways were associated with ≥1.5-fold DE genes found to significantly increase in response to 20 hr anoxia (Fig. 3A). A number of these pathways were related to (1) metabolism, (2) DNA damage and repair, as well as (3) PPAR signaling. The most interesting finding from our RNAseq analysis was that enzyme-encoding genes involved in glucose catabolism, including ‘citrate cycle TCA cycle’ and ‘pentose phosphate pathway’, were found to be enriched in our distinct up-regulated dataset at the transcript level. Details of this and other selected pathways are discussed below in the context of the expected physiological changes and metabolic demands placed upon the liver of anoxic turtles. Opposed to KEGG pathways that were found to significantly increase in response to anoxia, 50 KEGG pathways were associated with genes that were ≥1.5-fold DE down-regulated (Fig. 3B). Among these GS analysis enriched KEGG pathways, a number of pathways were related to (1) growth factor and cytokine signaling and (2) adhesion-related processes. Individual gene regulation patterns of select anoxia-responsive upregulated KEGG pathways are visualized in either DNA damage and repair (Fig. 4) or metabolism-related processes (Fig. 5). Both nucleotide excision repair (Fig. 4A) and mismatch repair (Fig. 4B) are displayed as examples of DNA damage and repair processes that are found to be upregulated in response to anoxia. The pentose phosphate pathway (Fig. 5A) and the citric cycle (Fig. 5B) have been shown as key metabolism-related anoxia-responsive pathways.

A visual presentation of anoxia-regulated DNA damage signaling and repair KEGG pathways.

Figure 4: A visual presentation of anoxia-regulated DNA damage signaling and repair KEGG pathways.

Pathways are shown for (A) nucleotide excision repair and (B) mismatch repair found to be up-regulated from distinct directional GS analysis. The fold change (logarithmic transformed) for the relevant genes are also plotted with green and red representing increase and decrease, respectively.
A visual presentation of anoxia-regulated DNA damage signaling and repair KEGG pathways.

Figure 5: A visual presentation of anoxia-regulated DNA damage signaling and repair KEGG pathways.

Pathways are shown for (A) pentose phosphate pathway and (B) citrate cycle KEGG pathways found to be up-regulated from distinct directional GS analysis. The fold change (logarithmic transformed) for the relevant genes are also plotted with green and red representing increase and decrease, respectively.


It is known that the complete deprivation of oxygen can present several significant, and often fatal, consequences to the cell. Besides the aforementioned stress on metabolism, oxygen deprivation is characterized by an altered redox environment, as well as decreased cellular pH resulting from lactic acid production. In the present study, we used RNAseq to explore the transcriptomic response of T. scripta elegans liver to periods of anoxia. As previously discussed, T. scripta elegans utilizes a number of different molecular approaches to adapt to a low oxygen environment and is able to avoid anoxia-induced cellular damage through the reprioritization and depression of metabolic processes, such that cellular energy demand can be met by anaerobic metabolism (Storey, 2007). While changes in gene expression have been previously observed in a targeted manner (Krivoruchko & Storey, 2015; Storey, 2007), in this study we attempt to discuss the transcriptome-wide changes that occur in response to anoxia in turtle liver tissue to elucidate the molecular pathways and functions that aid in surviving energy limited systems.

Reflective of the metabolically depressed state of the turtle in response to low oxygen, we found that the majority of genes were found to either decrease in expression or be unchanged from normoxia levels. Of the total 3,105 genes that showed statistically significant changes in gene expression between the two animal groups, 64 genes increased greater than 5-fold when compared to control values (Table 2), while 417 genes decreased to less than 20% of control values (Table S1). Although we did not explore the function of the 2,624 genes that did not fit our criteria to define genes as anoxia-response, we further investigated the 481 anoxia-responsive genes to determine whether anoxic-responsive changes in the transcriptome are reflective of known physiological responses.

Our gene set enrichment approach to identifying pathways that are enriched in anoxia-responsive genes (Fig. 3) revealed a number of pathways for future in-depth analysis. Interestingly, among the top 20 up-regulated pathways, 25% are involved in DNA repair and nucleotide metabolism (nucleotide excision repair, RNA polymerase, mismatch repair, RNA degradation and DNA replication), 25% are metabolism-related (maturity onset diabetes, adipocytokine signaling, glyoxylate and dicarboxylate metabolism, citric acid TCA cycle, and fatty acid metabolism), and 15% protein translation/degradation (ribosome, spliceosome and ubiquitin-mediated proteolysis) (Fig. 3A). Other pathways include circadian rhythm and immunodeficiency. Given the ATP-deprived cellular state that is imposed during anoxia, in a mechanism likely in response aimed towards energy conservation, an even greater number of pathways were found to decrease (Fig. 3B). Overall, our analysis indicates that T. scripta elegans may increase cellular capacity to (2) respond to and repair DNA damage, and (2) reprioritize anaerobic metabolism utilizing mechanisms as our top biological processes that may be dynamically regulated by the turtle to respond to anoxia and will be explored in greater detail in the sections below.

Oxidative stress and DNA damage repair

T. scipta elegans is able to survive several weeks of oxygen deprivation while overwintering by severely depressing their metabolic rate. However, upon recovery from anoxia, the re-introduction of oxygen brings about an overproduction of reactive oxygen species (ROS). The turtle is able to tolerate these, typically lethal, high levels of oxidative stress by maintaining an uncharacteristically high basal level of antioxidant defence mechanisms (Willmore & Storey, 1997; Krivoruchko & Storey, 2010).

In this study, the gene expression of nitric oxidase 1 (NOX1), a ROS generating enzyme, was determined to increase 5.63-fold following 20 hrs anoxia exposure (Table 2). The physiological function of NOX1 is not currently clear, but given its well-established oxidase activity, it likely involves the production of superoxide. Originally, NOX1 was described as an NADPH oxidase that stimulated mitogenesis (Arnold et al., 2001); however, subsequent studies collectively suggested NOX1 to be a driver of hypertension (Matsuno et al., 2005; Dikalova et al., 2005; Gavazzi et al., 2006). Interestingly, anoxia has been previously found to induce a ten-fold reduction in heart rate and blood flow (Overgaard, Gesser & Wang, 2007). This reduction in cardiac output better matches the requirements of a metabolically depressed state, while blood pressure is still largely maintained through an increase in vascular resistance (Overgaard, Gesser & Wang, 2007). It is possible that NOX1 may have a role in the vasoconstriction displayed in anoxic turtles. In addition to the turtle’s high level of antioxidant defense, our GS analysis suggests that enhanced mechanisms of DNA damage repair may also be an important component of anoxia survival (Fig. 4). Indeed, these repair mechanisms may be enhanced to either maintain genome stability throughout the anoxic exposure or to help protect DNA against the oxidative stress of reperfusion injury.

Critically, ROS is known to produce several forms of DNA damage, including modification of all nucleotide bases, deletions, frame-shifts, strand breaks, chromosomal rearrangements, and intra- and interstrand DNA cross-links (Gonzalez-Hunt, Wadhwa & Sanders, 2018). The ERCC1 (excision repair cross-complementation group 1) and XPF (i.e., ERCC excision repair 4, endonuclease catalytic subunit; a.k.a ERCC4) genes were both found to increase significantly in response to anoxia, among other genes that are involved in nucleotide excision repair, as per KEGG analysis (Fig. 4A; Table S1). Interestingly, together the ERCC1 and XPF genes encode the two subunits of the ERCC1–XPF nuclease. This complex plays an important role in repair of DNA damage and in maintaining genomic stability (Wang et al., 2011). The ERCC1–XPF complex is a core component of nucleotide excision repair and also facilitates inter-strand crosslink repair, pathways of double-strand break repair by homologous recombination and end-joining, as well as telomere length regulation (Wang et al., 2011). Interestingly, a previous study has shown that DNA repair-deficient ercc1 mutant mice displayed numerous features of accelerated aging and reduced lifespan (Vermeij et al., 2016). Furthermore, ercc1-deficient mice have been shown to have increased levels of 8-oxoguanine (i.e., the most common form of DNA lesion resulting from ROS) in their DNA (Hsia et al., 2003). Other central genes involved in DNA damage repair include the sliding clamp proliferating cell nuclear antigen (PCNA) that acts as an important platform for recruiting components of the DNA damage response and that was also found to significantly increase in response to anoxia, as suggested by both DE and GS analyses (Fig. 4B; Table S1) (Mailand, Gibbs-Seymour & Bekker-Jensen, 2013). Together, these findings suggest that DNA damage genes are actively up-regulated to repair oxidative-induced DNA damage, and potentially maintain genome integrity throughout oxidative stress. Indeed, the possibility of anoxia-induced increase in DNA damage repair capacity to help combat increased oxidative stress-related DNA damage has been previously proposed (Krivoruchko & Storey, 2010) and has been further supported by our findings.

Reprioritization of anaerobic metabolism

In aerobic animals, the complete catabolism of all fuel sources relies on oxidative phosphorylation, in which oxygen plays a crucial role as the final electron acceptor. When oxidative phosphorylation is inhibited, anaerobic pathways of energy production become the sole energy producers of the cell (Storey, 2007). However, simply increasing the rate of anaerobic glycolysis has several limitations, including poor ATP yield, toxic end products (lactic acid), and a quick depletion of fuel sources. Therefore, most aerobic life forms are fairly intolerant of oxygen limitation and sustain significant damage and/or death when subjected to anoxia.

In contrast, T. scripta elegans has developed several mechanisms to deal with the challenges posed by utilizing anaerobic glycolysis as their primary ATP producing pathway (Storey, 2007). Although increases in glycolytic rate have been reported for almost all organs of T. scripta elegans within the first five hrs of oxygen deprivation, there is a very rapid glycolytic inhibition in the liver occurring within the first hour of anoxia exposure. This reflects the need for glycogenolysis to be directed toward glucose export, supplying the fermentative fuel to be used by other organs (Kelly & Storey, 1988). As a result, a significant metabolic reorganization and reprioritization takes place in the liver in response to anoxia. Our RNAseq analysis has identified a number of differentially expressed metabolic enzymes and anoxia-responsive metabolic pathways. For example, significant alterations in the expression of a number of glycolytic enzymes in response to 20 hrs of anoxia exposure were suggested by the RNAseq DE results. This includes a 17-fold increase in the glucokinase binding protein, GKRP (Table 2); encoding a protein that modulates glucokinase (GK; activity and location of gluconeogenesis as glucose levels decline. In conditions of low glucose, GKRP relocates GK to the nucleus where it is inactive. Increasing amounts of glucose prompts the GKRP to rapidly release GK to return to the cytoplasm (Choi et al., 2013). Such a mechanism may help to dynamically control liver glucose metabolism in response to anoxia.

The transcript levels of glycolytic enzymes, phosphofructokinase (PFK; and aldolase (ALDO;, were also found to significantly increase in response to 20 hrs anoxia (Fig. 5A). Specifically, PFK significantly increased 2.77-fold, whereas ALDO increased 4.89-fold when compared to normoxic control values. Although changes in transcript abundance are not necessarily reflective of changes in enzyme function or activity, previous studies have reported significant increases in the activities of these two enzymes in response to anoxia in T. scripta elegans liver. One study documented significant changes in the enzyme activity of PFK upon prolonged exposure to anoxia, demonstrated a 3-fold decrease in the I50 value (i.e., the inhibitor concentration reducing enzyme velocity by 50%) for citrate, and found a 1.5-fold increase in the Km value (i.e., the substrate concentration producing half maximal enzyme activity) for ATP during anoxia (Brooks & Storey, 1989). Corroborating our RNAseq analysis, the protein expression of ALDO has been previously shown to significantly increase by 1.6-fold in response to 20 hrs anoxia in T. scripta elegans liver (Dawson, Biggar & Storey, 2013). In that study, the increase in fructose-1,6-bisphosphate (FBP) levels known to occur during anoxia were found to be accompanied by both an increase in endogenous ALDO protein expression, along with higher affinity for FBP in turtle liver. Another study also reported that the relative levels of FBP in liver doubled after only 1 hr of anoxia in T. scripta elegans (Kelly & Storey, 1988). When considering the presence of large liver glycogen stores and general increase in anaerobic glycolysis seen during anoxic exposure, the increase in FBP levels suggests the regulatory importance of aldolase in anoxic liver tissue.

Our RNAseq analysis also identified a significant increase in the gene expression of citrate cycle enzymes that support succinate production, which was enriched by the GS KEGG analysis (Fig. 5B). In this regard, it has been established that anaerobic glycolysis with lactate as the end product has two limitations: (1) a low yield of ATP per glucose catabolized, and (2) significant cellular acidification. Although the turtle has unique mechanisms to help deal with the latter, several alternative end products to lactate have appeared in anoxia-tolerant species that provide enhanced ATP yield per glucose catabolized. Collectively, these include ethanol, alanine, succinate, propionate, and acetate (Muller et al., 2012). For marine invertebrates, during prolonged anoxia, carbon from glycolysis is directed into the reactions of succinate synthesis. The energy yield of direct fermentation of glucose to succinate yields a molecule of GTP per molecule glucose in addition to the two molecules that result from lactate production. Although a mechanistic exploration into the enzymatic production of succinate as an anoxia end product of glucose catabolism in T. scripta elegans has not yet been completed, previous studies have explored alternative glycolytic end products in response to prolonged anoxia in C. picta bellii (Buck, 2000). In response to a 28-day anoxic dive at 5 °C, succinate, alanine, and lactate concentrations were all found to increase significantly in liver tissue. Relative to levels of lactate, this study found that levels of succinate and alanine increased by 2 and 0.9% in liver, 0.3 and 0.04% in blood, and 0.6 and 0.07% in heart, respectively. Together, these findings suggest a role for succinate production in turtle liver tissue as a possible glycolytic alternative to lactate. Indeed, our RNAseq analysis supports the hypothesis that succinate accumulation may be a component of anaerobic metabolism (Fig. 5B). For example, at the level of the transcriptome, our analysis found a significant increase in the anoxic expression of succinyl-CoA-synthetase (SCS), which was found to be elevated 1.62-fold in comparison to normoxia values. Critically, the SCS enzyme is responsive for the catalysis of the reversible reaction of succinyl-CoA to succinate and the associated production of GTP. Collectively, these findings support the possibility of glucose catabolism to succinate as a possible alternative fate to lactate, yielding greater energy yield but is likely a minor anaerobic end-product compared to lactate.


Although much is known about the physiological responses of anoxic exposure in the turtle, T. scripta elegans, studies are beginning to elucidate the molecular underpinnings that coordinate these adaptations. As research begins to utilize strategies of transcriptomic and proteomic-wide analysis to explore models of metabolic depression and extreme environmental adaptations, we will begin to better understand the molecular mechanisms that are in place to support survival. For example, this study highlighted a number of anoxia-responsive molecular pathways that are likely critical to navigating anoxia survival, such as DNA damage repair and metabolic reprogramming. Furthermore, the possibility of succinate production being utilized as an additional energy production pathway in addition to lactate metabolism was also highlighted as part of this study. Ultimately, this study has set the stage for future research to navigate and document the many insights that have been uncovered.

Supplemental Information

Enrichment boxplots

DOI: 10.7717/peerj.8144/supp-1

Raw RNAseq statistics

DOI: 10.7717/peerj.8144/supp-2

GO biological process enrichment statistics

DOI: 10.7717/peerj.8144/supp-3

GO molecular function enrichment statistics

DOI: 10.7717/peerj.8144/supp-4

KEGG pathway enrichment statistics

DOI: 10.7717/peerj.8144/supp-5