Development-specific transcriptomic profiling suggests new mechanisms for anoxic survival in the ventricle of overwintering turtles

For almost every animal, oxygen is required for life, but variation in survival time in the absence of oxygen, also called anoxia, exists both between and within species. The most anoxia-tolerant tetrapod is the North American pond turtle, Chrysemys picta, or the painted turtle. It can survive ∼170 days submerged in anoxic water at 3°C (Jackson et al., 2000; Odegard et al., 2018; Ultsch and Jackson, 1982), and for at least 30 h at 20°C (Johlin and Moreland, 1933). Painted turtles use this ability to withstand harsh winters in the northern regions of their geographical range, where they naturally experience severe hypoxia or anoxia within ice-covered ponds (Reese et al., 2004a; Ultsch, 1989). This extreme anoxia tolerance is achieved through metabolic suppression, utilization of large tissue glycogen stores, especially in the liver (Buck et al., 1993; Herbert and Jackson, 1985; Jackson, 1968, 2002; Warren and Jackson, 2007), and defense of body-fluid pH against lactic acidosis by utilizing bone as a buffer (Jackson and Ultsch, 1982; Warren and Jackson, 2008, 2017).

Although all organs in the turtle must tolerate anoxia, the heart is unique because it must simultaneously decrease metabolic rate to levels sustainable by anaerobic metabolism, and prevent circulatory arrest. Thus, a fine balance between cardiac functional suppression and maintained cardiac output must be achieved to survive anoxia. During anoxia, functional suppression is achieved by decreasing both heart rate and contractility, allowing continued function at minimal energy cost (Farrell et al., 1994; Hicks and Farrell, 2000; Hicks and Wang, 1998; Shi et al., 1999; Shi and Jackson, 1997; Wasser et al., 1990a,b). Even though functional shifts have been characterized, concordant transcriptomic changes during and following anoxia are poorly understood (Keenan et al., 2015).

As in most vertebrates, hypoxia/anoxia tolerance in the painted turtle is stage-dependent. Unlike adults, hatchling painted turtles can survive just 40 days when submerged in anoxic water at 3°C, which is similar to the anoxia tolerance of both hatchlings and adults of other pond turtle species (Dinkelacker et al., 2005b). This pattern of survival differs from that of other vertebrate species, of which neonatal animals show increased tolerance to hypoxia/anoxia compared with adults (Adolph, 1969). When viewed from an ecological perspective, the pattern observed in the painted turtle is not entirely surprising, because hatchlings overwinter in terrestrial nests where they avoid seasonal anoxia (Packard and Packard, 1993, 2001) and, instead, survive subzero temperatures (Churchill and Storey, 1992a,b; Dinkelacker et al., 2005b; Packard and Packard, 1993, 2001, 2004; Rubinsky et al., 1994; Storey et al., 1988). This developmental difference in anoxia tolerance provides a unique opportunity to compare transcriptomic shifts in response to anoxia from animals with a common genomic background.

The limited work characterizing the anoxia response in hatchling turtles has focused on lactic acid buffering, leading to the hypothesis that lower bone mineral content and, therefore, buffering capacity, prevent hatchlings from surviving anoxia for as long as adults (Packard and Packard, 2004; Reese et al., 2004b). However, to our knowledge, no studies have examined cardiac responses of hatchling turtles to anoxia. Buffering capacity may not be the only developmental advantage; adult-specific, anoxia-induced changes in cardiac function may prevent cardiac failure and allow the adult to survive. Furthermore, hatchlings exhibit higher lactate levels after 40 days of anoxia, reflecting a higher anaerobic metabolic rate than adults, which is consistent with the hypothesis of adult-specific, cardiac metabolic suppression in response to anoxia (Reese et al., 2004b).

We tested the hypothesis that transcriptomic patterns during anoxia and recovery in adult ventricular tissue would reveal upregulation of pro-anoxia survival pathways, while hatchling ventricular tissue would reflect active growth and development pathways that either must be maintained or cannot be downregulated. We quantified anoxia-induced mRNA expression using RNA sequencing (RNA-seq) in the ventricle of adult and hatchling painted turtles submerged in 3°C anoxic water for 20 days, followed by 5 days of reoxygenation. We found candidate genes potentially important for anoxic survival in adults and others that may be maladaptive for anoxic survival in hatchlings, possibly exemplifying a selective mechanism akin to antagonistic pleiotropy that favors the evolution or maintenance of extreme terrestrial overwintering behavior in natal nests.

In order to avoid removing breeding females from their natural populations, this study utilized only male painted turtles, Chrysemys picta (Schneider 1783). Because painted turtles exhibit temperature-dependent sex determination (Janzen and Morjan, 2002), this practice could not have impacted the genomic variation of the study (i.e. sex-linked genes do not exist in this species). Adult male painted turtles (25 males, 186–425 g) were captured (May–August 2015) and removed from six localities within 160 km of the greater St Louis area, MO, USA (permits: Illinois, NH15.5803; Missouri, 16516) (Table S1). Turtles were housed in fiberglass aquaria filled with dechlorinated, St Louis municipal water (18–22°C) with access to basking platforms bathed with 10 W UV and 60 W incandescent heating lights that followed a local Missouri photoperiod (updated weekly). They were fed ReptoMin^® three times a week ad libitum and chicken liver once per week.

Neonatal turtles were produced by collecting eggs laid by turtles injected with oxytocin (June–July 2015) in the laboratory (Tucker et al., 2007). All gravid painted turtles (N=9, range=354.5–529.6 g) were captured from four of the same or nearby localities where the adult males were taken. Thirty-four eggs (4–13 eggs from each location) were incubated (64–76 days) at 25°C in moistened vermiculite to produce all male offspring (Janzen and Morjan, 2002; Schwarzkopf and Brooks, 1985). These hatchlings (N=34, mean±s.e.m. mass=4.91±0.12 g, range=3.70–6.71 g) were transferred to Styrofoam containers filled halfway with autoclaved sand and held at 20°C without feeding until used in the study. Hatchlings were misted weekly with autoclaved, deionized water. All housing and animal procedures were approved by the Saint Louis University Institutional Animal Care and Use Committee protocol 2198.

Hatchling turtles were introduced to water and acclimated to 20°C by first placing them in a small, darkened, plastic aquarium (22.86×15.24×16.51 cm), the base of which had been removed and replaced with plastic grating. This smaller aquarium was then suspended within a larger, darkened, glass aquarium (30.48×121.92×45.72 cm, ≈170 liters) that later housed the adults during experimentation. This arrangement allowed the adults and hatchlings to be physically separated, but still share the same water. The aquarium system was partially filled with circulating, aerated, St Louis municipal water. Ten days after the hatchlings were introduced to water, adult turtles were added to the outer glass aquarium and the water level was adjusted so the body mass to accessible water ratio was the same for both developmental stages (51.8 g l⁻¹ H₂O). After four more days of acclimation at 20°C, the water temperature was decreased daily by 2°C and, upon reaching 10°C, water temperature was decreased by 1°C daily until reaching 3±0.1°C. A YSI model 72 proportional temperature controller was used in concert with an immersion water heater and a circulating water pump to maintain the water temperature. The cold acclimation, anoxia exposure and tissue sampling were all conducted in a 3°C walk-in environmental room.

Turtles were distributed into treatment groups to deliberately retain homogeneity of genetic variance. Adults were distributed with four to five different populations represented at each time point and no population was represented more than once. Hatchlings were distributed so that five or more populations were represented at each time point and no more than two clutch-mates were represented in the same treatment. Adult and hatchling turtles were sampled at three time points: after 47 days in 3°C water (3°C control; N=5 adults, 8 hatchlings), after 20 days of anoxic submergence at 3°C (anoxia; N=5 adults, 8 hatchlings) and after 5 days of recovery at 3°C (50% recovery; N=4 adults, 10 hatchlings). There were no differences in the masses of turtles sampled between time points within each developmental stage (adults: N=4–5, one-way ANOVA, P=0.88; hatchlings: N=8–10, one-way ANOVA, P=0.88).

After the control turtles were sampled, the water level was increased and the turtles for the anoxia treatment were submerged underneath plastic grating placed approximately 5–10 cm below the water's surface, preventing access to the air or gas space above. The water was then bubbled with nitrogen gas to displace the dissolved oxygen. The top of the tank was covered with glass and sealed with silicone sealant to prevent gas exchange between the atmosphere and the aquarium. The glass lid included a 1 inch (2.54 cm) hole and was plugged with a cored, rubber stopper. All wires associated with controlling and monitoring experimental conditions were threaded through the stopper. Nitrogen bubbling was continuous during the anoxic period and was monitored (D200 DO meter, YSI, Yellow Springs, OH, USA) throughout the 20 days of anoxia (0.01–0.05 ppm). Nitrogen gas was allowed to escape through small holes in the rubber stopper. The anoxic turtles were sampled through a door fashioned into the plastic grating without allowing them to breathe air.

After experiencing anoxia, the remaining turtles were transferred to a smaller aquarium filled with recirculating, aerated water that maintained the same biomass to water ratio (51.8 g l⁻¹ H₂O) at 3°C. The turtles had free access to air during recovery. To monitor plasma lactate levels during recovery, tail vein blood samples (0.05–0.10 ml) were taken daily from only the adult turtles. Such monitoring was impossible in the hatchlings because of their small size. When plasma lactate levels decreased to ∼50% of those during anoxia after 121–126 h of reoxygenation, all remaining turtles, including hatchlings, were sampled.

Turtles were sampled immediately after removal from the aquarium at 3°C. When anoxic turtles were sampled, the neck of each turtle was clamped underwater to prevent oxygenation of the blood, followed immediately by decapitation and pithing of the brain in air. During sampling of hatchling turtles, the plastron was removed using surgical scissors and the ventricle was immediately excised with sterilized spring scissors, briefly blotted on sterile surgical gauze to remove residual blood, and flash-frozen on liquid nitrogen-cooled freeze clamps. The plastron of the adult turtle was removed using a bone saw, and the ventricle was immediately excised with sterile surgical scissors and trisected on sterile aluminium foil. Each ventricular section was briefly blotted on a sterile surgical sponge to remove residual blood and flash-frozen. All samples were stored at −80°C prior to RNA extractions.

Blood samples from both adults and hatchlings were obtained in duplicate or triplicate by draining blood that accumulated in the pericardial cavity either into a heparinized syringe (adults) or into heparinized microhematocrit tubes (hatchlings). A subsample of adult blood was transferred to a microhematocrit tube and spun with hatchling blood samples in a microhematocrit centrifuge (IEC Model MB IM-173, Damon, Needham Heights, MA, USA) for 3 min. The hematocrit was read, the tubes were scored and snapped, and the plasma was recovered for immediate lactate and glucose measurements (YSI 2300 Stat Plus). The remaining adult blood was centrifuged at 9600 g (AccuSpin Micro 17, Fisher Scientific, Germany) and the plasma was recovered. All remaining plasma was transferred to microcentrifuge tubes and flash-frozen in a slurry of dry ice and ethanol.

Frozen ventricular muscle (32.7±5.9 mg, mean±s.e.m.) was mechanically homogenized to a fine powder using sterile and RNase-free zirconium mortars and pestles (Cryogrinder Kit, OPS Diagnostics, Lebanon, NJ, USA) that were precooled with liquid nitrogen. To maintain consistency in powdering methodology between developmental stages, a randomly-selected fragment of frozen ventricle was powdered from the adults (range=15.9–78.7 mg), while the whole hatchling ventricle (7.1–19.6 mg) was powdered. The powder was then transferred to sterile, RNase-free cryovials precooled with dry ice. Room temperature TRIzol reagent (20 µl mg⁻¹) (Life Technologies, Burlington ON, Canada) was added to the powder and immediately vortexed until it was completely suspended. Extractions were subsequently performed according to the manufacturer's guidelines. The resulting RNA pellets were resuspended in 1 mol l⁻¹ sodium citrate, pH 6.4±0.2 (The RNA Storage Solution, Life Technologies, Carlsbad, CA, USA), passed through a Zymo Clean and Concentrator-5 kit with DNase I (Zymo Research, Irvine, CA, USA), and eluted into 1 mol l⁻¹ sodium citrate. Final RNA concentrations were measured on a Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and purity (RNA integrity number mean=8.8, range=7.9–9.3) was confirmed using an Agilent 2100 bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).

cDNA library construction and mRNA sequencing were carried out on RNA samples from four adults and four hatchlings per treatment and time point at the McDonnell Genome Institute at Washington University (St Louis, MO, USA). The TruSeq stranded mRNA kit (Illumina, San Diego, CA, USA) was used to prepare polyadenylated-selected RNA-seq libraries of paired-end 2×150 bp reads (insert size≈200 bp). qPCR was used to amplify and normalize libraries prior to sequencing on an Illumina HiSeq 4000 platform. Samples were randomly distributed across five lanes of the flow cell, and there were no significant differences in sequencing depth (mean=97.6 million reads, range=75.4–116.3 million reads) between treatment (two-way ANOVA, P=0.33) or developmental (P=0.20) groups. Raw reads are publicly available on NCBI (project PRJNA526071).

TruSeq LT adapter contamination at 3′ ends was trimmed away using Cutadapt software (v1.4.2) and trimmed reads longer than 75 bp were retained and aligned to the Chrysemys picta bellii RefSeq Genome assembly sequence (v3.0.3) using Tophat 2 (v2.1.1) with guidance from the C. p. bellii reference annotation (v3.0.3) (Kim et al., 2013; Martin, 2011). Alignments were then sorted and indexed using Samtools (v1.3.1) for visualization against the genome using IGV software (v2.4.9) (Li et al., 2009; Robinson et al., 2011).

The Cufflinks 2 suite of tools (v2.2.1) was used to determine differential expression of annotated genes and transcripts. Briefly, Cufflinks 2 was used to assemble transcriptomes by generating normalized FPKM values (fragments per kilobase of transcript per million mapped reads) of previously annotated genes for each sample (Trapnell et al., 2012). Using Cuffmerge 2, sample transcriptomes were merged with the reference annotation to create a master reference transcriptome. Cuffdiff 2 calculated geometrically normalized mean FPKM values for every gene at each time point and developmental stage, and generated differential gene expression comparisons using a linear model assuming a normal/Gaussian distribution (Roberts et al., 2011). Fold-change comparisons were performed by log₂ transforming the quotient of mean FPKMs, adjusting P-values for multiple test comparisons using the Benjamini–Hochberg correction (Benjamini and Hochberg, 1995). Genes were excluded from differential expression analysis if the sum of all FPKM values across all collection time points for a given gene was less than 50. Significant differential expression was defined by an adjusted P-value ≤0.05, and a fold-change ≥1 or ≤−1 between mean FPKMs. Statistical comparisons of gene groups were analyzed and figures were generated using R v3.4.0 (https://www.r-project.org/) and the following packages: VennDiagram, ggplot2 and ComplexHeatmap (Chen, 2017; Gu et al., 2016; Wickham, 2009). Gene ontology (GO) analysis was carried out using the PANTHER classification system (Mi et al., 2017).

A two-factor ANOVA was used to analyze the effect of treatment and developmental stage on plasma lactate, glucose and hematocrit, and post hoc analysis was completed using Tukey HSD pairwise tests. Lactate and glucose values were square root-transformed to maintain homogeneity of variance. Differences were considered significant when P-values were <0.05.

After 20 days of anoxia, plasma lactate increased considerably from 2.6±1.0 to 48.3±6.7 mmol l⁻¹ (mean±s.e.m., pairwise t-test, P<0.001) in adults, and from 0.6±0.1 to 57.3±2.9 mmol l⁻¹ (P<0.001) in hatchlings (Fig. 1A). After 5 days of recovery, plasma lactate levels at both developmental stages decreased substantially by 54.4% and 43.1% to 22.0 and 32.6 mmol l⁻¹ in adults (P<0.001) and hatchlings (P<0.001), respectively, but remained elevated compared with controls in both groups (adult: P<0.001, hatchling: P<0.001). Plasma lactate did not differ between adults and hatchlings at any time point (control: P=0.39, anoxia: P=0.50, recovery: P=0.13). One adult and five hatchling turtles were non-responsive after the 20 days of anoxia and were declared dead.

Plasma glucose shifted ontogenetically (Fig. 1B). In adults, plasma glucose levels (Fig. 1B) were similar to control levels (2.93 mmol l⁻¹) after 20 days of anoxia (P=0.91) and after 5 days of recovery (P=0.74), whereas in hatchlings glucose increased from 1.67 to 11.0 mmol l⁻¹ after 20 days of anoxia (P<0.001) and continued to rise after 5 days of recovery, reaching 17.2 mmol l⁻¹ (P<0.001).

Blood hematocrit levels were affected by treatment (P=0.03). In both hatchlings and adults, hematocrit averages trended higher after 20 days of anoxia (Fig. 1C), although differences were not significant (adult: P=0.95, N=5; hatchling: P=0.73, N=8). At each time point, the adult mean was higher than the hatchling mean; however, the difference between the two stages was only significant during recovery (P=0.03; adult: N=4, hatchling: N=10).

The FPKM sum of 50 was considered the optimum threshold where simultaneously genes with low expression (mean FPKM <2.08) were filtered out, but those that showed significant differential expression remained. Genes were ranked by overall transcript abundance, which was determined by ranking the FPKM sum (N=12; 4 transcriptomes at three sampling time points) for each developmental stage. The top 13 most abundantly transcribed genes under all conditions were protein-coding regions of the mitochondrial genome in both developmental stages. The expression levels of these genes did not vary during anoxia or recovery or across developmental stage (median FPKM range=8032.0–47,306.9).

Of the most abundantly transcribed protein-coding genes from the nuclear genome, 91 ranked in the 100 most abundant at both developmental stages, with the remaining ranking in the top 200. Of the top 109 genes from adults and hatchlings combined, 71 were related to translational regulation (Table S2), where 65 were ribosomal protein (R-protein) genes and six were involved in translational control (UBA52, EEF1A1, EEF1G, EEF2, FAU, CRIBP-like). Of the remaining 38 genes, 13 were associated with the contractile apparatus and excitation–contraction coupling (ACTA1, ACTB, ACTG1, ACTN2, DES, MYL3, MYL7, MYL10, MYH15, TNNT2, TNNI3, TPM4, TPM1) and five were related to the mitochondrial electron transport chain and ATP synthesis (GAPDH, LDHB, NDUFA4, NME-like, SLC25A4). Three were associated with iron metabolism (FTL-like, FTH, FTH-like), three encoded globins (HBAA, HBB, MB) and three were associated with protein turnover (GNB2L1, UBB, UBC). The remaining 11 genes were related to a variety of other biological processes and do not easily group according to their function.

Analysis of ranked transcript abundance detected distinct development-specific patterns in R-protein gene expression in response to anoxia (N=65). Not surprisingly, because of their high abundance (median=1794.0 FPKM, range=176.2–6178.9 FPKM), the 76 individual R-protein genes did not reach the significance threshold in response to treatment, which requires a log₂ change of 1; the total R-protein abundance accounted for 15.4% of total FPKMs. However, further analysis revealed treatment-dependent expression patterns (Fig. 2). A Kolmogorov–Smirnov test detected a difference between adult and hatchling distributions of R-protein transcript FPKM values (P<0.001, N=76). Therefore, these distributions were analyzed separately using Kruskal–Wallis nonparametric one-way ANOVAs. The distributions of R-protein transcript FPKMs in adults and hatchlings showed a treatment effect (P<0.001 and P=0.036, respectively, N=76). Dunn post hoc tests with a Benjamini–Hochberg correction revealed that 20 days of anoxia decreased FPKM values for R-protein genes in adults (P=0.002) and increased them in hatchlings (P=0.044). The FPKM values returned to control levels after 5 days of recovery in both adults (P=0.24) and hatchlings (P=0.91).

Mitochondrial R-protein gene expression showed far less abundance than nuclear R-protein gene expression (median FPKM=28.1, FPKM range=1.0–209.7). Although a Kolmogorov–Smirnov test detected a developmental effect (P<0.001) on the expression of mitochondrial R-protein genes, Kruskal–Wallis one-way ANOVAs did not detect an effect of treatment in either adults (P=0.76) or hatchlings (P=0.14).

Multidimensional scaling analysis of Euclidian distances for the catalog of genes showed that transcriptomes clustered by developmental stage (Goodness of fit=0.29, 0.29) (Fig. 3A). Genes were analyzed individually to determine how many expression changes were attributable to development versus treatment. Out of the catalog of 11,072 genes, a total of 1260 (11.4%) experienced major changes in transcript abundance (Fig. 3B). Of these, developmental stage affected the expression of 1175 genes. Prior to anoxia, 567 genes differed in expression between adults and hatchlings. After 20 days of anoxia, 754 differed, of which 362 were different prior to anoxia. After 5 days of recovery, 703 genes differed between adults and hatchlings, of which 427 were also different prior to anoxia. Overall, at 3°C, 346 genes differed significantly between adult and hatchling turtles during all treatments.

Anoxia and reoxygenation accounted for substantial changes in expression of 237 genes, of which only 31 changed in both hatchlings and adults (Fig. 4). In adults, 130 changed uniquely, with 76 changing uniquely in hatchlings (Fig. 4C). After 20 days of anoxia, only three genes greatly increased expression levels in both adults and hatchlings and only four genes did so during recovery (Fig. 4D,F). Over the same period, seven genes decreased expression, whereas six decreased expression after 5 days of recovery at both developmental stages (Fig. 4E,G). The genes BTG2 and SFRP4 were downregulated during anoxia and recovery in both adults and hatchlings, indicating that only 18 out of 237 genes changed similarly owing to treatment across developmental stage.

After 20 days of anoxia, adults exhibited 64 uniquely upregulated genes and 41 uniquely downregulated genes, yet hatchlings only exhibited half or fewer uniquely upregulated and downregulated genes (Fig. 4D,E). A different pattern occurred during recovery: in adults, 39 unique genes increased expression and 45 decreased expression, whereas, for hatchlings, 43 unique genes were upregulated and 23 were downregulated (Fig. 4F,G). Overall, hatchlings had a larger proportion of increased expression; 67 of 107 genes were upregulated, with only 40 being downregulated during either anoxia or recovery.

Of the 237 protein-coding genes with differential expression during anoxia or reoxygenation, 48 were upregulated in adults to levels higher than observed in hatchlings (Table S4). Twenty-five genes were downregulated during anoxia in adults to levels lower than observed in hatchlings, while 33 were upregulated in hatchlings to levels higher than observed in adults (Table S5). Also, 14 genes were constitutively higher in transcript abundance in adults and were later upregulated in hatchlings during either anoxia or recovery. Twelve genes were always expressed at lower levels in adults and were later downregulated in hatchlings as a result of treatment (Table S6).

The present study is the first comparative analysis of anoxia-induced responses of ventricular transcriptomes from an animal that exhibits development-dependent anoxia tolerance. Adult and hatchling painted turtles showed some gene expression commonalities; few genes differed between the most abundant transcripts from both developmental stages. However, life stage remained the most important determinant of transcriptomic gene expression changes. Few genes changed significantly owing to anoxia or recovery, most of which were impacted by an interaction with developmental stage. Anoxia- and recovery-induced gene expression changes unique to developmental stage suggest candidate genes essential for surviving anoxia. Furthermore, developmental stage impacted R-protein gene expression patterns, suggesting that translational regulation is essential for anoxia tolerance in the adult turtle heart.

Lactate accumulation has been well characterized during anoxia in adult turtles and is often used as a proxy for anaerobic metabolic rate (Jackson and Heisler, 1982, 1983; Jackson et al., 2000; Warren and Jackson, 2007, 2017). Although ontogeny did not induce a statistically detectable effect (P=0.19, N=4–8; Fig. 1A), plasma lactate levels in hatchlings trended higher after 20 days of anoxia. Previous work showed that 40 days of anoxia increased hatchling lactate levels higher than adults (Reese et al., 2004b). Thus, the rate of anoxia-induced lactate accumulation is probably higher in hatchling turtles, but is probably too variable to detect with our small sample sizes so early in the submergence period.

This is the first study to report that hatchling painted turtles have greater circulating glucose concentrations than adults after 20 days of anoxia and 5 days of recovery (Fig. 1B), and may result either from decreased glucose utilization or, more likely, greater liver glycogenolysis compared with adults. Hatchling turtles utilize glycogen and mobilize glucose from the liver, as whole-body glycogen levels decrease after 40 days of anoxic submergence and plasma glucose increases during hypoxia (Costanzo et al., 2001; Dinkelacker et al., 2005a; Reese et al., 2004b). The hatchlings also show decreased liver glycogen content with increased plasma glucose during terrestrial freezing (Packard and Packard, 2004, 2005). In contrast, glucose levels in adults did not differ after 20 days of anoxia or 5 days of recovery in the present study, indicating development-dependent differences in glucose homeostasis during and following anoxia.

A comparison of transcriptomes at 3°C reveals a common dependence on genes with the highest transcript abundance. The highest-ranking protein-coding genes were represented in high abundance at both developmental stages, indicating that the most abundant protein-coding genes at this temperature in painted turtles do not change ontogenetically.

A portion of the most abundant genes from both developmental stages in painted turtles is also highly abundant in human ventricular tissue, suggesting that the mRNA transcriptome for this tissue is evolutionarily conserved. A comparison of painted turtle with human ventricular tissue detected common expression patterns of 41 out of the 89 (46%) most abundant protein-coding genes in the human ventricle (Table S3). Eight of 89 were mitochondrial protein-coding genes, suggesting that mitochondrial genes are highly expressed in ventricular tissue regardless of species (Melé et al., 2015). Of the remaining 33 commonly transcribed protein-coding genes from both turtles (adults and hatchlings) and humans, 12 were related to translation regulation (Fig. S1), 11 of which were ribosomal protein-coding genes (RPL8, RPL10, RPL19, RPL26L1, RPL27, RPLP1, RPS11, RPS12, RPS16, RPS18, RPS27A). High transcript abundance of ribosomal proteins is not surprising, as they are essential for basic translational regulation and are commonly expressed in all tissue types across vertebrates (Glisovic et al., 2008; Hsiao et al., 2001).

One of the most important findings concerns the patterns of R-protein gene expression after anoxic exposure, which suggest a mechanism for anoxia-induced translational suppression unique to adult turtle ventricles (Hochachka and Lutz, 2001; Keenan et al., 2015) (Fig. 2). In adult turtles, overall R-protein gene expression declined after 20 days of anoxia at 3°C (P<0.001, N=76), but showed no difference from the control time point after 5 days of recovery (P=0.24), indicating restoration to pre-anoxic levels. The decrease in overall abundance of R-protein transcripts suggests a decrease in abundance of proteins required for translation, and therefore may constitute a mechanism for suppression of both translation and metabolism in adult turtle hearts. Anoxia decreases protein synthesis in the turtle ventricle (Bailey and Driedzic, 1996), which could derive from arrested translation owing to rRNA depletion with a decline in total RNA content. A previous study of anoxic painted turtles at 19°C showed a decrease in total RNA in the ventricle during anoxia (Keenan et al., 2015). However, such differences were not found in the present study. Instead, the present study suggests that at 3°C, R-protein levels decrease during anoxia in adults, which would be a cause of anoxia-induced translational arrest. This mechanism for translational arrest would not be unique to turtles; stress-induced suppression of R-protein gene expression has been observed in stress-resistant yeast, which shows a rapid decrease in R-protein transcription levels in response to environmental stresses including hydrogen peroxide-induced oxidative stress and heat shock (Gasch et al., 2000; Warner, 1999).

In stark contrast, hatchling turtles experienced an anoxia-induced increase in R-protein transcript abundance (P=0.015) and a return to control levels after 5 days of recovery (P=0.91), suggesting no modulation of translational activity in the ventricle during anoxic exposure as a means of metabolic suppression. Rather, the increase in R-protein expression may indicate an increase in translational activity in response to anoxia, which would be maladaptive under this energetically stressful condition (Fig. 2). Further investigation into R-protein abundance is needed to confirm that changes in R-protein gene expression influence protein levels during anoxia.

The most important factor affecting gene expression was developmental stage. Initial analysis of the gene catalog (N=11,072) showed clear grouping of all transcriptomes by developmental stage, indicating that adult and hatchling turtles exhibited constitutive differences in transcriptomes (Fig. 3A). Furthermore, analysis of individual genes showed that over the course of the three time points, the expression of 1175 genes significantly differed between the two developmental stages, while only 237 genes showed expression changes owing to treatment (Fig. 3B). The importance of developmental stage in gene expression is not a novel finding, as ontogenetic gene expression changes have been the focus of many studies (Gellon and McGinnis, 1998; Riggs and Podrabsky, 2017; Rougvie, 2001); however, this finding supports the argument that whole transcriptome comparisons across developmental stages can be used as a discovery tool to characterize changes in tissue development.

Developmental stage not only characterized fundamentally different transcriptomes, but also interacted with treatment effects, suggesting developmental-specific responses to treatments reflect developmental-specific anoxia tolerance. Of 237 differentially expressed genes caused by treatment (Fig. 4D–G), only 18 changed similarly in adults and hatchlings during anoxia or recovery. Therefore, adults and hatchlings not only exhibit transcriptomic differences at 3°C, but also different ventricular transcriptomic responses to anoxia and reoxygenation. Hence, protein expression and, consequently, ventricular function also may differ between adults and hatchlings before and during anoxia. Turtle ventricles appear to rely less on anoxia-induced gene expression changes per se, but, instead, develop their anoxia tolerance as they age by increasing the expression of constitutively adaptive genes and decreasing the expression of maladaptive ones.

By assuming the patterns of adult gene expression reflect an anoxia-tolerant, cardioprotective phenotype, and hatchling gene expression patterns reflect an anoxia-intolerant one, we have utilized a comparative approach to identify gene candidates that potentially characterize survival. These genes fall into three categories: (1) genes that are uniquely upregulated in adults; (2) genes that are inherently expressed at higher levels in adults, and are induced by anoxia in hatchlings (these genes are considered constitutively adaptive or protective); and (3) genes that are downregulated in adults, but remain elevated in hatchlings (these genes are considered maladaptive or injurious).

A total of 93 genes were upregulated in adult turtles after anoxia or recovery and may play a protective role during anoxic survival (Table S4). Of these genes, 42 were uniquely upregulated from control values in adults to levels significantly higher than observed in hatchlings after 20 days of anoxia. Eight of these genes were downregulated in hatchlings, further indicating a development-specific response to anoxia and reoxygenation (Table 1). Of these, Ficolin 2-like (FCN2-like), a signaling molecule involved in the lectin complement pathway, showed a 5-fold increase in expression in adults and a 2-fold decrease in hatchlings. Although Ficolin 1 and 3 are expressed in leukocytes and monocytes, respectively, FCN2 function is best understood in human liver, where it is secreted into the plasma to help initiate the lectin complement pathway (Kilpatrick and Chalmers, 2012). LYZ, a gene involved in lysozyme activity, TUBA8, a gene coding for a cytoskeleton protein involved in GTPase activity, and PATE3, an uncharacterized protein, exhibited this same pattern of expression. LYZ gene expression further increased after 5 days of recovery, as also seen in the hypoxia-tolerant freshwater fish Megalobrama amblycephala, which experiences increased LYZ gene expression during hypoxia and reoxygenation (Chen et al., 2017). The unique upregulation of these genes may allow adult turtles to survive anoxic conditions.

Of the 42 uniquely upregulated genes in adult turtles, two were involved in vascular function, indicating that anoxia induced changes in not only ventricular tissue, but probably also ventricular vasculature (Table 1). Both ACTA2, an aortic actin isoform, and CSRP2, a gene involved in smooth muscle development, increased expression during anoxia to higher levels in adults than in hatchlings. The differing expression levels of these genes could also reflect developmental differences in ventricular vascularization, which, to our knowledge, has not been studied.

After 5 days of recovery, 26 genes were uniquely upregulated from control values in adults compared with hatchlings (Table S4) and are, therefore, viewed as promoting survival during recovery in the anoxia-tolerant phenotype. Of interest is S100A1, which encodes for a Ca²⁺-binding protein that plays an important role in Ca²⁺ handling during cardiac function (Duarte-Costa et al., 2014; Wright et al., 2009) and may be protective upon reperfusion. During recovery, S100A1 increased expression from a mean FPKM of 53.0 to 128.9 in adults but remained at low expression levels in hatchlings (mean FPKM=34.71) (Table 1). Functionally, S100A1 modifies excitation–contraction coupling by affecting sarcolemmal Ca²⁺ flux through indirect modification of L-type Ca²⁺ channels and Na⁺/Ca²⁺ exchanger activity (Most et al., 2005; Reppel et al., 2005), and by affecting Ca²⁺-induced Ca²⁺ release (CICR) from the sarcoplasmic reticulum (SR) through direct modification of ryanodine receptors and SR calcium ATPase (SERCA2a) (Prosser et al., 2008; Remppis et al., 2002; Schaub and Heizmann, 2008; Wright et al., 2008). In painted turtles, the role of SR Ca²⁺ and S100A1 in ventricular CICR remains unclear. Previous work in turtle and trout ventricular myocytes suggested a lack of SR Ca²⁺ involvement in resting Ca²⁺ transients or twitch force (Cros et al., 2014; Galli et al., 2006a,b), but that simulated anoxia increased diastolic Ca²⁺ (Wasser and Heisler, 1997). An increase in SERCA2 activity during recovery could decrease diastolic Ca²⁺ in adult ventricular myocytes, which may allow larger Ca²⁺ transients and, therefore, a restoration of ventricular contractility in the face of acidosis-induced decreases in contractility (Fanter et al., 2017). Further work is needed to demonstrate that recovery-induced S100A1 gene expression alters protein expression and that turtle ventricular Ca²⁺ handling is modified.

Genes that increased abundance in both developmental stages after either 20 days of anoxia or 5 days of recovery could defend cardiac function during anoxia and reperfusion. Only 16 genes were commonly upregulated in response to treatment, three of which increased during anoxia (Table S4). One example is UCP3, which encodes uncoupling protein 3 (Table 1), a mitochondrial inner membrane protein that affects mitochondrial reactive oxygen species (ROS) production in skeletal muscle (Nabben et al., 2008; Vidal-Puig et al., 2000). After 20 days of anoxia, UCP3 expression increased 5.6-fold in adults and 2.3-fold in hatchlings. Increased UCP3 mRNA expression during anoxia may reflect a protective mechanism where UCP3 translation is initiated immediately upon reperfusion and, therefore, attenuates ROS production.

Based on comparative analyses, 20 genes may be constitutively adaptive or protective in adults. These genes all remained unchanged in adults, but were upregulated in hatchlings during anoxia or recovery to levels similar to those in adults (Table S6). Seven of these genes became upregulated in hatchlings after 20 days of anoxia, three of which coded for the histone proteins Histone H1, Histone H1.11L-like and Histone H2B 8 (Table 2). Increased expression of histone-encoding genes suggests that hatchlings have anoxia-induced changes in DNA structure and organization that, in turn, might change gene expression (Fan et al., 2005). Because expression of histone genes in adults did not change during anoxia, and their expression levels are inherently more abundant than in hatchlings, it is possible that elevated histone H1 and H2B 8 gene expression is characteristic of an anoxia-tolerant transcriptome. To our knowledge, this is the first study to implicate constitutively elevated expression of histone mRNA as a characteristic of anoxia tolerance in the ventricle, reflecting an advantage of the comparative approach we used.

Transcripts significantly downregulated in both adults and hatchlings after 20 days of anoxia or 5 days of recovery point to genes that may be maladaptive or injurious during anoxia. Six genes were commonly downregulated in response to anoxia, five in response to recovery and two during both treatments (Table S5). One of these 13 genes, RNA Transcription, Translation and Transport Factor (RTRAF; previously known as C14orf166), is downregulated 4.14-fold during anoxia in adults and 2.79-fold in hatchlings (Table 3). This protein could affect formation of mature tRNA, arguing for another possible mechanism of translational regulation during anoxia exposure (Popow et al., 2014). Interestingly, this gene remains downregulated during recovery in hatchlings, but returns to control levels in adults, which suggests another ontogenetic difference in translational regulation.

Another downregulated gene was ANKRD1, which encodes the cardiac ankyrin repeat protein (CARP), a transcription factor that interacts with sarcomere proteins such as desmin and titin (Table 2). Upon α-adrenergic stimulation, CARP is relocalized to the nucleus, where it is suspected to modulate gene expression (Boriek and Mohamed, 2012; Maeda et al., 2002; Miller et al., 2003; Zhong et al., 2015). ANKRD1 silencing also disrupts sarcomere integrity and attenuated hypertrophy-induced NPPA mRNA expression (Chen et al., 2012; Zhong et al., 2015), indicating it might affect both contractile function and ANP signaling during anoxia.

After 20 days of anoxia, 18 genes were downregulated in adults to levels substantially lower than in hatchlings. Genes in this category could be important for growth and development in hatchlings, but maladaptive for anoxic survival, leading to their suppression in adults (Table S5). One example, NPPA-like (Table 3), was downregulated in adults (from 107.4 to 16.7 FPKM), but upregulated in hatchlings (from 46.9 to 134.7 FPKM). NPPA-like is similar to human atrial natriuretic peptide precursor A, which is one of 12 proteins with the highest level of enriched expression in the proteome of the human heart (Melé et al., 2015). In mammals, ANP is secreted from the myocardium to induce vasodilation of blood vessels and renal excretion of both Na⁺ and water in response to high blood pressure and volume (Song et al., 2015). During anoxia, adult painted turtles decrease renal glomerular filtration rate, likely an important component of overall metabolic suppression (Warburton and Jackson, 1995). The adult-specific reduction in NPPA-like gene expression during anoxia may prevent ANP secretion, further limiting tubular secretion of Na⁺ and water, and contribute to renal metabolic suppression. In contrast, the significant increase in NPPA-like gene expression in hatchlings might be maladaptive for anoxic survival because it could increase renal metabolic rate.

Cysteine-rich angiogenic inducer 61 (CYR61), also known as CCN1, was downregulated during anoxia in adults, yet remained elevated in hatchlings (Table 3). CYR61 expression is important during cardiac development in mice (Mo and Lau, 2006). Under control conditions, hatchlings expressed CYR61 at levels 5.23-fold higher than adults, which showed a 3.1-fold and a 2.5-fold decrease in expression during anoxia and recovery, respectively. Hatchlings showed no change in CYR61 expression, but always had much higher expression levels than adults. We hypothesize that CYR61 is important for cardiac development in hatchlings and must either remain elevated during anoxia or cannot be downregulated. In adults, downregulating this gene may simply reflect anoxia-tolerant transcriptional suppression of a non-essential gene as a way of promoting cardiac survival.

Several genes involved in mediating oxidative stress were also uniquely downregulated in adults during anoxia, including two GSTM1-like genes that code for glutathione S-transferase (Table 3). Although control FPKM values were similar across developmental stage, only adults had major decreases in GSTM1-like mRNA expression during anoxia. Glutathione S-transferase (GST) aids in preventing oxidative stress-induced apoptosis of cardiomyocytes by helping detoxify ROS (McBride et al., 2005; Röth et al., 2011). Both the turtle heart and brain produce less anoxia-induced ROS than mammals, while reoxygenation-induced ROS levels in turtle neurons do not appear to differ from those in normoxic cells (Bundgaard et al., 2018; Milton et al., 2007; Pamenter et al., 2007). Furthermore, after 24 h of recovery, turtle heart glutathione levels were not oxidized, an indication of little ROS formation (Willmore and Storey, 1997). Adult turtles may suppress ROS production during reoxygenation and, therefore, have little need for GST-mediated ROS handling and GSTM1-like activity. In contrast, ROS may be produced in hearts of hatchlings and may pose a greater threat to hatchling survival than adults. It is also possible that hatchlings simply cannot downregulate these genes.

Genes upregulated during anoxia and recovery in hatchlings, but remaining at lower levels in adults, all potentially contribute to the inability of hatchlings to survive protracted anoxia and may even be considered maladaptive for survival. After 20 days of anoxia, 42 genes were uniquely upregulated in hatchling turtles, of which 30 increased to levels higher than observed in adults. Of these 30 genes, five encoded for different hemoglobin isoforms – HBAD, HBAA, HBE1-like and two HBB transcripts – four of which remained upregulated during recovery (Table 3). Additionally, PPOX, ALAS2 and ALAS2-like were upregulated during anoxia in hatchlings, and all but ALAS2 remained so during recovery. All eight of these genes play important roles in heme biosynthesis; the ALAS2 enzyme controls the rate-limiting step in heme biosynthesis and PPOX controls the final step (Sawicki et al., 2015). In mammals, ALAS2 levels are upregulated during hypoxia, supporting hypoxia-induced heme biosynthesis (Hofer et al., 2003). Elevated levels of heme have been observed in infarcted mice hearts and ALAS2 over-expression induces both increased heme content and associated oxidative stress in murine cardiomyoblasts (Sansbury et al., 2014; Sawicki et al., 2015). This suggests that hatchlings may somehow be more vulnerable to oxidative stress through upregulation of heme biosynthesis during and post anoxia, which may contribute to decreased survival during reoxygenation. In stark contrast, the adults showed no change in expression of any of these eight genes during both anoxia and recovery, maintaining levels considerably lower than those observed in hatchlings. It is important to note, also, that we have presumed these mRNA changes were occurring in the muscle, but we cannot rule out that they occurred in erythrocytes. Although mature erythrocytes typically arrest RNA and protein synthesis, developing red blood cells from both newts and chickens synthesize hemoglobin (Cameron and Prescott, 1963; Grasso et al., 1977).

Hatchlings also upregulated carbonic anhydrase 1 (CA1) gene expression during anoxia, which could have implications for cardiac pH regulation in vivo. Carbonic anhydrases catalyze the reversible hydration of CO₂ to carbonic acid, which dissociates into H⁺ and HCO₃⁻ under physiological pH, thereby regulating total CO₂ and H⁺ levels (Maren, 1967). Carbonic anhydrases have been well studied in the mammalian heart and are believed to selectively facilitate Na-coupled HCO₃⁻ transport as a mechanism for intracellular pH (pH_i) regulation (Villafuerte et al., 2014). If the levels of CA1 protein track those of the mRNA, then the capacity for HCO₃⁻ flux might also be upregulated in anoxic hatchlings. Such a response would help the cell defend intracellular pH only if the CA1 were externally localized on the sarcolemma, by enhancing HCO₃⁻ influx; however, if the CA1 were localized internally, it would acidify the sarcoplasm by enhancing HCO₃⁻ efflux and likely decrease contractility. Interestingly, both snapping turtle embryos and yearlings exposed to hypercapnia defend cardiac pH_i while extracellular pH decreases (Shartau et al., 2016). A specific role of CA1 in this response has not been elucidated in these species. For comparison, adult painted turtles do not defend pH_i during hypercapnic acidosis under normoxic or anoxic conditions (Jackson et al., 1991). Cardiac pH_i regulation has not been characterized in hatchling painted turtles.

From an evolutionary perspective, it is puzzling that painted turtles would maintain large numbers of stage-specific, potentially maladaptive genes, such as those downregulated in the adults during or following anoxia, while upregulated or maintained in the hatchlings (see above). These include the genes encoding the R-proteins (discussed earlier), some of the most highly expressed genes in the turtle ventricle. However, considering that these potentially maladaptive genes are rarely, if ever, under negative selection for a hatchling, the picture becomes clearer. Hatchling turtles typically overwinter in terrestrial nests and can face subzero temperatures for months, rather than enduring anoxic conditions underwater that adults experience (Churchill and Storey, 1992a,b; Dinkelacker et al., 2005b; Packard and Packard, 1993, 2001, 2004; Rubinsky et al., 1994; Storey et al., 1988). This implies that fitness costs associated with anoxia-induced upregulation/downregulation for hatchlings must be severe and that such costs are absent in the adults. These genes can be viewed as antagonistically pleiotropic because they confer fitness benefits to terrestrially overwintering hatchlings yet would have adverse fitness consequences if expressed in aquatically overwintering adults. Consequently, we suggest an innovative explanation for resolving this long-standing conundrum: avoiding antagonistic pleiotropy across life stages has favored the evolution, or at least maintenance, of the extreme terrestrial overwintering behavior of neonatal painted turtles.

Transcriptomic profiling of the painted turtle ventricle showed that ontogeny plays an important role in anoxia-induced changes in gene expression. A total of 1175 genes differed significantly between developmental stages across all time points. In contrast, after 20 days of anoxia and 5 days of recovery, only 237 genes changed, many of which were unique to each developmental stage. Adult painted turtles can survive anoxia four times longer than hatchlings, and this resilience may result, in part, from these observed developmental differences in ventricular gene expression. Our results revealed candidate genes that may characterize adult-specific anoxia tolerance and implicated many others as potentially maladaptive for hatchling survival during anoxia. Such patterns, which are akin to antagonistic pleiotropy across developmental stages, suggest a possible novel selective framework supporting the evolution or maintenance of terrestrial overwintering by hatchlings to avoid the substantive fitness costs of expressing those genes during anoxic conditions. Overall, these findings demonstrate the manifold power of applying a comparative approach combining two developmental transcriptomes with the same genomic background.

We thank Dr Claire Riggs for her useful discussion, and Craig Hill for his sampling assistance.