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

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 N₂ 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/.

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 log₁₀(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 log₂ 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.

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.

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. S1–S8).

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.

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; 2.7.1.2) 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; 2.7.1.11) and aldolase (ALDO; 4.1.2.13), 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 I₅₀ value (i.e., the inhibitor concentration reducing enzyme velocity by 50%) for citrate, and found a 1.5-fold increase in the K_m 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.