Developmental programming of DNA methylation and gene expression patterns is associated with extreme cardiovascular tolerance to anoxia in the common snapping turtle

Abstract

Background

Environmental fluctuation during embryonic and fetal development can permanently alter an organism’s morphology, physiology, and behaviour. This phenomenon, known as developmental plasticity, is particularly relevant to reptiles that develop in subterranean nests with variable oxygen tensions. Previous work has shown hypoxia permanently alters the cardiovascular system of snapping turtles and may improve cardiac anoxia tolerance later in life. The mechanisms driving this process are unknown but may involve epigenetic regulation of gene expression via DNA methylation. To test this hypothesis, we assessed in situ cardiac performance during 2 h of acute anoxia in juvenile turtles previously exposed to normoxia (21% oxygen) or hypoxia (10% oxygen) during embryogenesis. Next, we analysed DNA methylation and gene expression patterns in turtles from the same cohorts using whole genome bisulfite sequencing, which represents the first high-resolution investigation of DNA methylation patterns in any reptilian species.

Results

Genome-wide correlations between CpG and CpG island methylation and gene expression patterns in the snapping turtle were consistent with patterns observed in mammals. As hypothesized, developmental hypoxia increased juvenile turtle cardiac anoxia tolerance and programmed DNA methylation and gene expression patterns. Programmed differences in expression of genes such as SCN5A may account for differences in heart rate, while genes such as TNNT2 and TPM3 may underlie differences in calcium sensitivity and contractility of cardiomyocytes and cardiac inotropy. Finally, we identified putative transcription factor-binding sites in promoters and in differentially methylated CpG islands that suggest a model linking programming of DNA methylation during embryogenesis to differential gene expression and cardiovascular physiology later in life. Binding sites for hypoxia inducible factors (HIF1A, ARNT, and EPAS1) and key transcription factors activated by MAPK and BMP signaling (RREB1 and SMAD4) are implicated.

Conclusions

Our data strongly suggests that DNA methylation plays a conserved role in the regulation of gene expression in reptiles. We also show that embryonic hypoxia programs DNA methylation and gene expression patterns and that these changes are associated with enhanced cardiac anoxia tolerance later in life. Programming of cardiac anoxia tolerance has major ecological implications for snapping turtles, because these animals regularly exploit anoxic environments throughout their lifespan.

Introduction

The environment that an organism experiences in early life can have profound and long-lasting effects on their phenotype. This phenomenon, termed developmental plasticity, allows animals to permanently alter their morphology, physiology and behaviour in response to environmental signals [1]. In many cases, developmental plasticity provides organisms with a powerful mechanism to cope with environmental heterogeneity later in life [2]. However, unexpected or severe environmental stress during development can produce maladaptive phenotypes that increase disease susceptibility [3]. Despite the profound ecological implications of developmental plasticity, the underlying cellular and molecular mechanisms remain poorly defined.

Due to the profound health implications, most studies investigating developmental plasticity have focused on mammalian models of disease [4]. However, environmental variation during development is much more common in ectothermic animals, particularly oviparous species [5, 6]. These animals typically develop with little or no parental care and are routinely subjected to wide variations in abiotic factors such as temperature, water availability and atmospheric gases [7]. In particular, oviparous reptile nests can become severely hypoxic due to a progressive decline in nest oxygen tension from embryonic metabolism and microbial activity [8, 9]. The extent of hypoxia is nest-specific, but field estimates suggest reptilian eggs located farthest from the surface can be subjected to oxygen tensions as low as 11%, while those at the top of the nest remain at atmospheric oxygen (21%) [10]. Similar to other vertebrates, developmental hypoxia significantly alters turtle morphology and physiology, particularly at the level of the cardiovascular system [11,12,13,14,15]. Embryonic turtles exposed to hypoxia have different intrinsic heart rates and variable expression of receptors involved in cardiac regulation [11, 13, 16,17,18,19]. Furthermore, the effects of developmental hypoxia extend into juvenile and adult life, affecting cardiac performance and physiological traits [14, 15]. Of particular note, our recent study suggests juvenile turtles from hypoxic incubations possess cardiomyocyte specialisations that improve anoxia tolerance [20]. The programming of cardiac anoxia tolerance has major ecological implications for turtles, because many freshwater species, including Chrysemys picta, Trachemys scripta, and Chelydra serpentina, regularly engage in breath-hold dives that last several hours at warm temperatures, and they overwinter in anoxia for up to 5 months in ice-covered lakes [21, 22]. Even when metabolic rate and body temperature are taken into account, these freshwater turtles can survive anoxia 1000 times longer than a similarly sized mammal [23]. The maintenance of cardiac function is crucial for anoxia survival to ensure the delivery of nutrients and the removal of waste [24]. Therefore, early exposure to hypoxia may prime turtle heart physiology for a future life in anoxic environments.

The molecular mechanisms underlying cardiac programming in turtles are completely unknown but may involve epigenetic regulation of gene expression. Post-translational histone modifications and DNA methylation are the primary epigenetic marks shown to play a role in development and differentiation [25,26,27]. These marks regulate gene expression patterns, cell-fate decisions, and cellular physiology by altering DNA accessibility and chromatin structure. For example, trimethylation of histone H3 on lysine 4 (H3K4me3) at promoters is associated with gene activation, while trimethylation of lysine 27 on histone H3 (H3K27me3) is a repressive mark [28]. At least 70 different histone marks have been identified, each having unique effects on gene expression. The complexity of the histone code contrasts with the relative simplicity of DNA methylation, which is associated with transcriptional repression. DNA methylation is thought to inhibit transcription by interfering with transcription factor (TF) binding, though TF binding might reciprocally inhibit DNA methylation [29, 30]. Moreover, histone modifications and DNA methylation are interdependent, so de novo DNA methylation patterns laid down during embryogenesis help set the stage for maintenance of DNA methylation patterns and histone modifications through repeated cell divisions and into postnatal life [31, 32].

DNA methylation is a particularly stable, long-term mark that might be subject to environmental modification during development [33]. The most common mark is methylation of cytosines adjacent to guanines (i.e., CpG dinucleotides). Individual CpGs are typically methylated, while CpGs in clusters, called CpG islands (CGIs), are usually, though not always, found in an unmethylated state. The impact of CpG and CGI methylation on gene expression also depends upon their location within the genome. Recent work, for instance, has shown that enhancers and silencers display different patterns of CpG methylation and that orphan CGIs can act as potent enhancers [34,35,36]. This is on top of the classical observation that 60–70% of promoters contain CGIs [37].

Developmental hypoxia is known to alter DNA methylation and gene expression patterns in mammals, and the molecular signature is associated with cardiac abnormalities in adulthood [38, 39]. Therefore, programming of cardiac anoxia tolerance in snapping turtles may be achieved by similar mechanisms. Very little is currently known about DNA methylation landscapes in reptiles, because prior studies have almost exclusively measured global DNA methylation levels. We found one study that examined spatial patterns using MeDIP-Seq in the painted turtle, Chrysemys picta [40]. Key observations were that CpG distribution is bimodal in turtle promoters, as in other vertebrates [41], and that there is differential CpG methylation between hatchling ovaries and testes, including methylation differences in putative sex-determining genes. While MeDIP-Seq provides an overview of the methylation landscape at an affordable cost, it is an enrichment-based technique with shortcomings in terms of quantitatively measuring methylation levels and presenting a biased representation of the genome [42]. More importantly, we could not find a single study describing the most fundamental relationships between DNA methylation and gene expression patterns in reptiles.

In this study, we hypothesised that developmental hypoxia alters DNA methylation and gene expression patterns in turtles and that these patterns are associated with greater cardiac anoxia tolerance later in life. Snapping turtles take 9 to 18 years to reach sexual maturity, which makes it impractical to study developmental programming in adults. Instead, we tested for effects that persist in juvenile turtles months after their embryonic exposure to hypoxic conditions. To directly test these hypotheses, we first assessed cardiac performance during 2 h of acute anoxia in juvenile turtles previously exposed to normoxia (21% oxygen: N21) or hypoxia (10% oxygen: H10) during embryonic development. Next, we measured DNA methylation patterns in heart ventricles from the same cohorts using whole genome bisulfite sequencing (WGBS), the “gold standard” for DNA methylation analyses, as well as gene expression patterns using RNA-Seq. These experiments represent the first high-resolution investigation of DNA methylation patterns in any reptilian species. As hypothesized, developmental hypoxia increased juvenile turtle cardiac anoxia tolerance and programmed CpG and CGI methylation and gene expression patterns. DNA methylation and gene expression were broadly correlated at a genome-wide scale (e.g., genes with higher methylation at their promoters displayed lower expression, while those with lower promoter methylation displayed higher expression). In addition, genes that were differentially methylated between turtles from normoxic and hypoxic incubations were significantly more likely to be differentially expressed. The results suggest developmental hypoxia can programme turtle cardiovascular phenotype, spanning from molecular to physiological levels, which has important ecological implications for species that exploit anoxic environments.

Results

Developmental hypoxia improves cardiac anoxia tolerance

Body and heart masses of juvenile turtles used for in situ studies of cardiovascular physiology are provided in Table 1. Acute exposure to anoxia caused a progressive bradycardia (i.e., decreased heart rate) in both experimental groups (Fig. 1A), but the magnitude of this response was significantly greater in N21 (34 ± 6%) vs. H10 (20 ± 10%) turtles. A decrease in heart rate in low oxygen environments is a key feature of the “diving reflex”, which aids in the conservation of oxygen stores in air breathing vertebrates. In the N21 group, bradycardia was associated with a progressive reduction in systemic blood flow ($\dot{Q}_{{{\text{Sys}}}}$) and pulmonary blood flow ($\dot{Q}_{{{\text{Pul}}}}$) (Fig. 1B, C), while systemic stroke volume ($V_{{{\text{S}},{\text{Sys}}}}$) and pulmonary stroke volume ($V_{{{\text{S}},{\text{Pul}}}}$) remained relatively constant (Fig. 1E, F). The reduction in pulmonary blood flow ($\dot{Q}_{{{\text{Pul}}}}$) in N21 turtles during anoxia was proportionately greater than the reduction in systemic blood flow ($\dot{Q}_{{{\text{Sys}}}}$), leading to an increase in the right-to-left (R–L) shunt of blood from the pulmonary to the systemic circulation (Fig. 1H). Turtles are able to physiologically control the outflow of blood through the pulmonary artery vs. systemic arteries (i.e., left and right aortas), because they have a three chambered heart with a single ventricle that is only partially divided by vertical and horizontal septa. An increase in R–L shunting recirculates systemic venous blood and bypasses the pulmonary circuit, while greater left-to-right (L–R) shunting recirculates blood through the pulmonary circuit. Changes in shunting may allow more efficient regulation of blood gases during periods of activity vs. rest [43]. Despite a significant reduction in total blood flow ($\dot{Q}_{{{\text{Tot}}}}$) in N21 turtles (Fig. 1D), there was only a small non-significant reduction in cardiac power output (Fig. 1J), while mean ventricular pressure remained relatively constant (Fig. 1I).

Table 1 Body and heart masses of juvenile snapping turtles exposed to normoxia (N21) or hypoxia (H10) during embryonic development

Fig. 1

Effects of acute anoxia and reoxygenation on haemodynamic variables from N21 and H10 turtles. Turtles from the N21 (red circles, n = 6) and H10 (blue squares, n = 5) cohorts were subjected to 120 min of anoxia followed by 30 min reoxygenation. A Heart rate ($f_{{\text{H}}}$), B systemic blood flow ($\dot{Q}_{{{\text{Sys}}}}$), C pulmonary blood flow ($\dot{Q}_{{{\text{Pul}}}}$), D total blood flow ($\dot{Q}_{{{\text{Tot}}}}$), E systemic stroke volume ($V_{{{\text{S}},{\text{Sys}}}}$), F pulmonary stroke volume ($V_{{{\text{S}},{\text{Pul}}}}$), G shunt distribution ($\dot{Q}_{{{\text{Shunt}}}}$), H shunt ratio (${{\dot{Q}_{{{\text{Pul}}}} } \mathord{\left/ {\vphantom {{\dot{Q}_{{{\text{Pul}}}} } {\dot{Q}_{{{\text{Sys}}}} }}} \right. \kern-\nulldelimiterspace} {\dot{Q}_{{{\text{Sys}}}} }}$), I mean ventricular pressure ($P_{{{\text{Vent}}}}$), and J cardiac power output. Values are mean ± SEM, asterisks (*) indicate statistically significance difference between N21 and H10 groups, dollar ($) and psi (Ψ) symbols denote a significant difference between that data point and pre-anoxic levels (time zero) in the N21 and H10 groups, respectively (p ≤ 0.05)

Apart from the “diving reflex” (i.e., bradycardia), other cardiovascular responses in H10 turtles were quite distinct from N21 turtles. Surprisingly, the anoxic bradycardia in H10 turtles was not associated with any changes in systemic ($\dot{Q}_{{{\text{Sys}}}}$) or pulmonary ($\dot{Q}_{{{\text{Pul}}}}$) blood flow or the R–L shunt, which all changed in N21 turtles. This meant that systemic ($V_{{{\text{S}},{\text{Sys}}}}$) and pulmonary ($V_{{{\text{S}},{\text{Pul}}}}$) stroke volumes were significantly elevated in H10 turtles during acute anoxia (Fig. 2). As a result of the elevated stroke volume, mean ventricular pressure and cardiac power output was maintained during 2 h of anoxia in H10 turtles (Fig. 1I, J). Therefore, the H10 group maintained higher blood flows, systemic stroke volume, and heart rate ($\dot{Q}_{{{\text{Sys}}}}$, $\dot{Q}_{{{\text{Pul}}}}$, $V_{{{\text{S}},{\text{Sys}}}}$, $f_{{\text{H}}}$)and cardiac power output than the N21 cohort throughout the anoxic period (Figs. 1 and 2). In the N21 group, all haemodynamic variables reverted to normoxic levels after 30 min of reoxygenation (Fig. 1). For the H10 group, mean ventricular pressure was slightly depressed at the end of reoxygenation (Fig. 1I), and $V_{{{\text{S}},{\text{Pul}}}}$ remained elevated (Fig. 1F), while all other haemodynamic variables returned to normoxic levels (Fig. 1).

Fig. 2

Original traces of the effects of anoxia and reoxygenation on cardiac haemodynamic variables. Ventricular pressure ($P_{{{\text{Vent}}}}$), left aortic arch blood flow ($\dot{Q}_{{{\text{LAo}}}}$), left pulmonary artery blood flow ($\dot{Q}_{{{\text{LPa}}}}$) and heart rate ($f_{{\text{H}}}$) were measured in N21 (red lines) and H10 (blue lines) turtles during 10-min normoxia, 120-min anoxia, and 20-min reoxygenation

In addition to influencing responses to anoxia and reoxygenation, developmental hypoxia altered resting cardiovascular variables in snapping turtles, similar to previous reports [14]. While all haemodynamic variables fell within previously published in situ and in vivo values from Chelydra, Chrysemys, and Trachemys [14, 44], anaesthetised H10 turtles had significantly greater resting systemic blood flow (Q_Sys)than N21 turtles, leading to a larger R–L shunt, elevated total blood flow (Q_Tot) and elevated cardiac power output (Fig. 1, pre-anoxic levels). All the other haemodynamic variables were similar between experimental groups.

Embryonic hypoxia programs transcriptome-wide patterns of gene expression

Transcriptome-wide patterns of gene expression were investigated in 7- and 9-month-old snapping turtles previously exposed to hypoxia (H10, n = 8) or normoxia (N21, n = 8) during embryonic development. Within the hypoxic cohort, turtles had two distinct cardiac phenotypes; normal-sized (n = 4) and enlarged (n = 4) hearts, relative to their body size. Gene expression within both cohorts was found to be significantly affected by age, relative heart size, and embryonic oxygen concentration. Firstly, oxygen concentration during embryogenesis altered expression of 151 genes in juvenile turtles: 75 genes were up-regulated and 76 genes were down-regulated in ventricles from the H10 group, relative to the N21 group (Table 2). Ninety-seven genes displayed significant oxygen concentration by age interactions (Table 3) and 13 of these genes were also influenced by the main effect of oxygen concentration. Finally, 256 genes were differentially expressed between ventricles from normal-sized vs. enlarged hearts (47 of these genes were among the genes listed above). A total of 131 genes were up-regulated in ventricles of enlarged hearts, while 125 genes were down-regulated (Table 4). Altogether, there were 443 differentially expressed genes.

Table 2 Genes that were differentially expressed between ventricles from juvenile snapping turtles exposed to normoxia (N21) or hypoxia (H10) during embryonic development

Table 3 Genes that displayed significant oxygen concentration by age interactions in ventricles from juvenile snapping turtles exposed to normoxia (N21) or hypoxia (H10) during embryonic development and sampled at 7 months or 9 months of age

Table 4 Genes that were differentially expressed between ventricles from juvenile snapping turtles that had normal-sized or enlarged hearts relative to their body size

Hierarchical clustering of these genes by expression pattern showed separation of normal-sized from enlarged hearts (i.e., the two deepest branches in the dendrogram in Fig. 3A). There was also separation between younger and older turtles (the next deepest branches in the dendrogram). Finally, two distinct clusters contained the N21 and H10 groups from 7-month-old turtles (separation of red and blue branches in top half of the dendrogram). Overall, this pattern of clustering reflected clear expression differences between N21 and H10 groups.

Fig. 3

Gene expression patterns in ventricles of 7-month-old and 9-month-old snapping turtles that had been incubated in normoxic or hypoxic conditions as embryos. Samples from the N21 cohort are shown with red circles and red lines, while the H10 cohort is shown with blue circles and blue lines. A Heatmap of RNA-Seq expression values for the 443 genes that were significantly affected by oxygen concentration during embryogenesis (Table 2), the oxygen concentration by age interaction (Table 3), and/or differed between ventricles from juvenile snapping turtles that had normal-sized or enlarged hearts relative to their body size (Table 4). Reverse transcription of total RNA and qPCR with rigorous standard curves were used to measure expression of 16 genes (B–Q) identified as differentially expressed in the RNA-Seq study. Significant differences between oxygen groups or oxygen by age interactions were confirmed for 14 genes (panels B through O), but not for 2 genes (panels P and Q). Expression levels are least squares means (± 1 SE) for each oxygen group at 7 and 9 months of age

Differentially expressed genes were enriched for several GO terms important for cardiac function and/or remodeling (Fig. 4). For GO Biological Processes, this included 8 differentially expressed genes that play a role in sarcomere organization, 33 genes that play a role in biological adhesion/cell adhesion, and 118 genes involved in signal transduction (Fig. 4; Table 5). For GO Cellular Components, 9 differentially expressed genes form collagen trimers, 11 genes are part of the Z-disc, 41 genes are found in the extracellular space, and 50 genes are part of the extracellular region (Fig. 4; Table 6). Several genes across different GO categories are candidates that might play a role in promoting cardiac anoxia tolerance in the H10 group.

Fig. 4

Gene Ontology terms significantly enriched among genes that were differentially expressed in ventricles from juvenile snapping turtles that had been incubated in normoxic or hypoxic conditions as embryos. Enrichment analysis was carried out on all 443 genes that were affected by oxygen concentration during embryogenesis (Table 2), the oxygen concentration by age interaction (Table 3), and/or those genes that differed between ventricles from turtles that had normal-sized vs. enlarged hearts relative to their body size (Table 4). The number of genes for each GO term is represented by the size of the circle, while the FDR adjusted p value is shown by the color of the circle

Table 5 GO biological process terms and genes that were enriched among differentially expressed genes in ventricles from juvenile snapping turtles

Table 6 GO Cellular Component terms and genes that were enriched among differentially expressed genes in ventricles from juvenile snapping turtles

We selected genes for qPCR validation from the GO categories described above based on their established role in influencing cardiac function and anoxia tolerance, including genes associated with heart defects in humans or other species, genes involved in calcium signaling or mitochondrial function, and/or genes that regulate expression of other genes. Overall, differential expression was confirmed for 14 of 16 genes examined (Table 7; Fig. 3B–Q). Some genes, such as DDIT4L and WNT11, were expressed at consistently lower levels in the H10 group compared to the N21 group at both ages (Fig. 3B, C). Other genes, such as ITGA11, MIPEP, MNAT1, PPIA, TNNT2, and TPM3, were reliably higher in the H10 group compared to the N21 group (Figs. 3D–I). Several genes displayed treatment by age interactions. For COL8A1, NCOA2, and SCN5A there was no difference at 7 months of age, but expression was higher in the H10 group than the N21 group at 9 months of age (Fig. 3J–L). For HIF1A and PFKFB1, there was no difference at 7 months of age, but expression was lower in the H10 group than the N21 group at 9 months of age (Fig. 3M, N). Another pattern was observed for HTRA3, which differed between treatment groups at 7 months of age, but not at 9 months of age (Fig. 3O). In contrast, CALR and SNTB1 did not differ between N21 and H10 groups at either age (Fig. 3P, Q).

Table 7 Results from a two-way ANCOVA for mRNA expression in ventricles from juvenile snapping turtles

Genome-wide correlation between DNA methylation and gene expression

DNA samples from ventricles of 9-month-old turtles were used for WGBS (n = 3 from N21 and n = 3 from H10). Given that DNA methylation landscapes have never been examined at a genome-wide scale in any reptile, basic patterns of DNA methylation were characterized before testing for differences between treatment groups. The draft snapping turtle genome contains approximately 142.4 million CpG dinucleotides. The distribution of CGIs in different genomic features was not random with respect to the proportion of the genome found in promoters, gene bodies, and intergenic regions: more CGIs were found in promoters (0 to − 1000 bp from the transcription start site) and in gene bodies than expected by chance, while fewer CGIs were found in intergenic regions (Table 8). For CpGs with sufficient read coverage (≥ 10 reads in 2/3 of replicates), levels of methylation were high (two thirds of CpGs across the genome were 75–100% methylated) and there were clear differences in methylation patterns among genomic features (Fig. 5). Intergenic regions (Fig. 5A) had a broader range of DNA methylation levels than did gene bodies (Fig. 5B). Intergenic regions had a higher proportion of CpGs with 0–75% methylation and a lower proportion of CpGs with 75–100% methylation than did gene bodies (Table 9). In other words, gene bodies were more heavily methylated than intergenic regions.

Table 8 Summary of the observed and the expected number of CGIs in different genomic features

Fig. 5

CpG methylation and gene expression patterns in the ventricles of juvenile snapping turtles. Probability density function for CpG methylation level in various genome features (A–F). A Intergenic regions (n = 12,840,311), B gene bodies (n = 3,426,421), C promoters (n = 68,168), D first exons (n = 40,226), E remaining exons (242,295), and F introns (n = 3,1561,179). Violin plots showing CpG methylation levels for various gene features as a function of gene expression level (G–K). Genes were divided into expression deciles with the lowest expressed genes in the first decile and the highest expressed genes in the tenth decile as shown on the X axis. The width of the violin corresponds to the probability density function for the methylation level shown on the Y axis. The bar in the middle of each violin shows the interquartile range and the white dot shows the median methylation level for each expression decile. Spearman’s rank correlations (ρ) between DNA methylation and gene expression were highly significant for all gene features (***p < 0.001). G Gene bodies (n = 3,426,421), H promoters (n = 68,168), I first exons (n = 40,226), J remaining exons (242,295), and K introns (n = 3,1561,179)

Table 9 Summary of the number of CpG sites by their genomic location and methylation level

There were also differences in CpG methylation patterns among gene features. Promoters (0 to − 1000 bp from the transcription start site or TSS) and first exons displayed a bimodal pattern of DNA methylation (Fig. 5C, D; Table 9), with a higher proportion of CpGs with 0–25% methylation (including unmethylated CpGs) and a lower proportion of CpGs with 75–100% methylation when compared to gene bodies, remaining exons, and introns (Fig. 5B, E, F; Table 9). That is, promoters and first exons displayed more variation in methylation levels and were less methylated on average than other exons and introns.

To test whether there was any relationship between CpG methylation and gene expression, genes expressed at a detectable level in turtle ventricles were divided into deciles based on expression levels with the first decile containing genes that displayed the lowest expression and the tenth decile containing genes that displayed the highest expression. There was a positive correlation between CpG methylation in gene bodies and expression levels (Fig. 5G). In contrast, CpG methylation in promoters and first exons was negatively correlated with gene expression (Fig. 5H, I). Remaining exons and introns displayed a positive correlation between CpG methylation and gene expression (Fig. 5J, K).

CpG methylation levels were plotted as a function of distance from TSSs to examine the methylation landscape of promoters at a finer spatial scale. Genes were divided into quintiles based on expression levels with the first quintile containing genes with the lowest expression and the fifth quintile containing genes with the highest expression. There was a clear negative correlation between methylation and gene expression levels (Fig. 6A). Genes in the first expression quintile exhibited slightly higher CpG methylation near the TSS vs. neighboring sites (i.e., a hill). In contrast, genes in the second through fifth expression quintiles exhibited progressively lower CpG methylation near the TSS (i.e., greater depth of the methylation valley with increasing expression). This valley spanned from roughly 1500 bp upstream to 1500 bp downstream of the TSS (Fig. 6A). A scatterplot of methylation levels for individual CpGs for genes in the fifth quintile showed a clear bimodal pattern centered on the TSS (i.e., most sites displaying 0% or 100% methylation) (Fig. 6B). This analysis demonstrated an inverse relationship between CpG methylation and gene expression: higher methylation at TSSs was associated with lower gene expression, while lower methylation at TSSs was associated with higher expression at a genome-wide scale.

Fig. 6

CpG methylation levels around TSSs as a function of gene expression level in ventricles of juvenile snapping turtles. A Genes were divided into expression quintiles with the lowest expressed genes in the first quintile and the highest expressed genes in the fifth quintile. Methylation level for each quintile was plotted as a function of the distance from the TSS. Lines were produced via LOESS curve fitting. B Scatter plot showing methylation level for each CpG as a function distance from the TSS for genes in the fifth quintile. The blue line is the same LOESS smoothed curve as shown in the previous panel. C CpG methylation level around TSSs for genes that were differentially expressed in ventricles of 9-month-old snapping turtles that had been incubated in normoxic (N21) or hypoxic (H10) conditions as embryos. Methylation levels are shown for genes that were up-regulated (red line) or down-regulated (green line) in hypoxic turtles compared to normoxic turtles

Embryonic hypoxia programs genome-wide patterns of CpG and CpG island methylation

Given that CpG methylation patterns were broadly correlated with gene expression in turtle hearts, fetal programming of DNA methylation could be driving hypoxia-induced differences in gene expression and physiology. The first step toward testing this hypothesis is to determine whether embryonic exposure to hypoxia caused differential DNA methylation. CpGs and CGIs were examined separately, because methylation patterns in mammals differ between isolated CpGs (heavily methylated) vs. CGIs (lightly methylated), as does the relationship of CpG and CGI methylation to gene expression.

Comparison of N21 and H10 groups revealed 74,016 differentially methylated CpGs out of 10,808,104 CpGs with sufficient coverage for analysis and difference > 25% and q < 0.01. Hypoxic conditions during embryogenesis induced hypermethylation of 38,428 CpGs and hypomethylation of 35,588 CpGs. Intergenic regions were not more or less likely to contain differentially methylated CpGs than expected by chance (Odds Ratio = 1.014, 95% CI = 0.999 to 1.029; Fisher’s Exact p = 0.066). However, differentially methylated CpGs were more likely to be found in promoters than expected by chance (Odds Ratio = 1.148, 95% CI = 1.057 to 1.245; Fisher’s Exact p = 0.001). In contrast, differentially methylated CpGs were less likely to be found in first exons (Odds Ratio = 0.707, 95% CI = 0.613 to 0.811; Fisher’s Exact p = 2e−07) or the remaining exons (Odds Ratio = 0.510, 95% CI = 0.481 to 0.541; Fisher’s Exact p = 2e−16).

Comparison of N21 and H10 groups revealed 6,666 differentially methylated CGIs (FDR < 0.05). Hypoxic conditions during embryogenesis induced hypermethylation of 3628 CGIs and hypomethylation of 3038 CGIs. Intergenic regions were more likely to contain differentially methylated CGIs than expected by chance (Odds Ratio = 1.136, 95% CI = 1.078 to 1.197; Fisher’s Exact p = 1.33e−6). In contrast, differentially methylated CGIs were less likely to be found in promoters than expected by chance (Odds Ratio = 0.738, 95% CI = 0.596 to 0.907; Fisher’s Exact p = 0.003). Although not statistically significant, a trend toward fewer differentially methylated CGIs was also observed in first exons (Odds Ratio = 0.815, 95% CI = 0.654 to 1.006; Fisher’s Exact p = 0.06) and the remaining exons (Odds Ratio = 0.847, 95% CI = 0.702 to 1.016; Fisher’s Exact p = 0.07).

Functional enrichment among differentially methylated genes

For GO analysis, differentially methylated genes (n = 1582) were defined as those containing ≥ 1 differentially methylated region (methylKit 200 bp sliding window with 50 bp steps) within their promoter (− 1000 bp from TSS) and/or gene body at a q < 0.001 (Additional file 1: Table S1). Differentially methylated genes were significantly enriched for numerous GO terms at a Bonferroni corrected p < 0.05 (Fig. 7). Eighteen terms were significant for GO Biological Process, including six GO terms that might be related to differences in the autonomic nervous system and bradycardia between N21 and H10 groups (Fig. 7; Additional file 2: Table S2). Among these, the highest level terms include “regulation of trans-synaptic signaling”, “regulation of nervous system development”, and “regulation of neuron differentiation” (Fig. 7; Additional file 2: Table S2). Thirty-one GO terms were significant for GO Cellular Component (Fig. 7; Additional file 2: Table S2). Several of these terms were also related to neuronal function, while other terms were related to cation channels that could play a role in positive ionotropic responses in the H10 group. Finally, GO Molecular Function contained 10 terms that complement GO Biological Process and Cellular Component terms (Fig. 7; Additional file 2: Table S2).

Fig. 7

Gene Ontology terms significantly enriched among genes that were differentially methylated in ventricles from juvenile snapping turtles that had been incubated in normoxic or hypoxic conditions as embryos. Enrichment analysis was carried out on 1582 genes that were affected by oxygen concentration during embryogenesis. The number of genes for each GO term is represented by the size of the circle, while the FDR adjusted p value is shown by the color of the circle

Correlation between hypoxia-induced DNA methylation and gene expression patterns

Having demonstrated that embryonic exposure to hypoxic conditions programmed differential methylation of CpGs and CGIs in juvenile turtle hearts, it was possible to test for relationships to hypoxia-induced differences in gene expression. Genes containing at least one differentially methylated region (as defined in the previous paragraph) were more likely to be differentially expressed than expected by chance (Odds Ratio = 1.558, 95% CI = 1.178 to 2.059; Fisher’s Exact p = 0.002). Genes that were both differentially methylated and differentially expressed between the N21 and H10 groups are listed in Table 10. Given the negative correlation between CpG methylation in promoters and gene expression (Fig. 5H) and the clear methylation signal centered on TSSs (Fig. 6A), finer scale CpG methylation patterns were examined for genes that were differentially expressed between the H10 and N21 groups at 9 months of age. Genes that were up-regulated and down-regulated by hypoxic incubation exhibited spatially distinct methylation patterns, particularly in the − 200 to − 1000 bp region of promoters (Fig. 6C). Although differences were not as stark as observed at a genome-wide scale, there was higher methylation in this region for genes that were downregulated by hypoxic incubation and lower methylation for upregulated genes. Taken together, these findings suggest that hypoxia-induced differences in CpG methylation in promoters and/or gene bodies contribute to hypoxia-induced differences in gene expression.

Table 10 Genes that were both differentially methylated and differentially expressed between ventricles from snapping turtles that were exposed to normoxia (N21) or hypoxia (H10) during embryonic development

Distal regulatory elements such as enhancers and silencers are another important factor driving gene expression patterns. Recent work in mammals has shown that orphan CGIs can act as enhancers. It is, therefore, possible that differential methylation of “CGI enhancers” could influence gene expression patterns in the snapping turtle. As a preliminary test of this hypothesis, we identified the closest gene to the 6666 differentially methylated CGIs described above. There were 4379 protein-coding genes near these sites. We then tested whether these genes were more or less likely to be differentially expressed between N21 and H10 groups. Genes closest to differentially methylated CGIs were less likely to be differentially expressed than expected by chance (Odds Ratio = 0.732, 95% CI = 0.578 to 0.927; Fisher’s Exact p = 0.005). Long-distance enhancer–promoter interactions that skip over the closest gene may explain this result (see discussion for in-depth consideration of this idea).

Putative cis-regulatory sequences for fetal programming of DNA methylation

HOMER2 was used to identify cis-regulatory elements that could be directing differential DNA methylation and/or differential gene expression between N21 and H10 groups. Proximal promoters (− 1000 bp to TSS) of the 443 differentially expressed genes were retrieved from the snapping turtle genome and analyzed for known and de novo sequence motifs. Known motifs for androgen response elements and glucocorticoid response elements were found in promoters at a higher frequency than expected by chance, but were not significant after FDR correction. Thirty-three (33) de novo sequence motifs were enriched in promoters (Additional file 3: Table S3). As expected for promoters, HOMER2 identified TATA boxes and downstream core elements (DCEs), which are recognized by TATA-binding protein and TFIID, respectively, at TSSs. The two highest scoring de novo motifs, which exceeded a stringent cutoff of 1e−10, were similar to sequences bound by transcription factors ZNF711 and RREB1. Several de novo sequence motifs are recognized by transcription factors (CEBPB, KLF10, MEIS1, RREB1, and RXRA) known to bind methylated DNA in mammals [29].

Given that CGIs can act as enhancers in other species, it is possible that differential methylation of CGIs could modulate their activity as enhancers in turtles. HOMER2 was, therefore, used to identify motifs within the 6,666 differentially methylated CGIs. A total of 42 known sequence motifs were enriched in differentially methylated CGIs with an FDR < 0.005 (Additional file 4: Table S4). Three transcription factors that bind to these motifs play a central role in mediating the effects of hypoxia on gene expression: i.e., HIF1A, ARNT (encodes Hif1-beta), and EPAS1 (encodes Hif2-alpha)-binding sites were all enriched in differentially methylated CGIs. Several transcription factors (CUX2, GABPA, GSC, PHOX2B, SMAD4, and ZEB2) that bind to enriched sequence motifs are associated with human syndromes that exhibit cardiovascular, mitochondrial, or autonomic nervous system defects (Additional file 4: Table S4). Several transcription factors (CRX, ZBTB33, SMAD4) are also known to bind methylated sites [29].

Discussion

The present study provides the first ever genome-wide analysis of DNA methylation patterns in a reptile and gives a detailed characterization of the methylation landscape across the snapping turtle genome. In addition, we show that developmental hypoxia is a potent environmental stimulus that alters snapping turtle DNA methylation patterns, which are associated with changes in gene expression and improved performance of the cardiovascular system during anoxia. Surprisingly, juvenile turtles from hypoxic incubations were able to maintain cardiac pumping capacity throughout 2 h of anoxia at 30 °C, which is a feat unsurpassed among vertebrates. To understand the molecular mechanisms underlying programmed differences in cardiac physiology, we carried out WGBS and RNA-Seq studies in ventricular tissue from the H10 and N21 cohorts. Developmental hypoxia programmed genome-wide methylation patterns at CpGs and CGIs. Furthermore, DNA methylation patterns were broadly correlated with gene expression at a genome-wide scale as well as for genes that were differentially expressed between normoxic and hypoxic turtles. Finally, we identified enriched DNA sequence motifs (i.e., putative TF-binding sites) in promoters of differentially expressed genes and in differentially methylated CGIs. By integrating this information, we develop a hypothetical model that links hypoxia during embryogenesis to persistent changes in DNA methylation and gene expression that are related to programmed differences in cardiomyocyte and cardiac physiology.

Effects of developmental hypoxia on the cardiovascular response to anoxia and reoxygenation

Similar to previous in situ [45] and in vivo [44] studies on turtles, acute anoxic exposure had negative chronotropic and inotropic effects in the N21 cohort. At the end of acute anoxia, systemic and pulmonary blood flow ($\dot{Q}_{{{\text{Sys}}}}$, and $\dot{Q}_{{{\text{Pul}}}}$) in N21 turtles was significantly reduced by 51% and 44%, respectively, leading to an increase in the R–L shunt. The reduction in total blood flow ($\dot{Q}_{{{\text{Total}}}}$) in the N21 cohort was achieved by a pronounced bradycardia, while stroke volume ($V_{{{\text{S}},{\text{Total}}}}$) remained unchanged, indicating a reduction in cardiac inotropy (contractility). These results align with previous work that demonstrates vagally mediated bradycardia in anoxic turtles at 21–25 °C [44, 46] as well as negative inotropy driven by intracellular acidosis, a reduction in calcium uptake, and energy depletion [20, 47]. Reduced cardiac activity during anoxia is characteristic of turtles and serves to lower ATP demand below the capacity for anaerobic ATP supply to restore energy balance [48]. Indeed, cardiac power output remained relatively stable during anoxia, indicating that ATP supply and demand were matched [48].

In contrast to the N21 cohort, H10 turtles maintained systemic and pulmonary blood flow ($\dot{Q}_{{{\text{Sys}}}}$, and $\dot{Q}_{{{\text{Pul}}}}$,) throughout the anoxic period. Maintenance of total blood flow ($\dot{Q}_{{{\text{Total}}}}$) was supported by a blunted bradycardia and an increase in stroke volumes ($V_{{{\text{S}},{\text{Sys}}}}$ and $V_{{{\text{S}},{\text{Pul}}}}$), while ventricular pressure and cardiac power output stayed constant. Surprisingly, in contrast to previous in vivo and in situ studies of anoxia-tolerant turtles (Trachemys scripta and Chrysemys picta) [44, 49, 50], 2 h of anoxia in H10 snapping turtles led to an increase in total blood flow ($\dot{Q}_{{{\text{Total}}}}$) and maintenance of cardiac pumping capacity. These findings are remarkable when considering the 30 °C body temperature of turtles in the current study compared to prior studies of animals held at 22 °C [44, 49, 50]. The only other vertebrate known to increase total blood flow ($\dot{Q}_{{{\text{Total}}}}$) and maintain pumping capacity during anoxia is the crucian carp, Carassius carassius, a species that can survive anoxia for months while remaining active [51]. Nevertheless, the crucian carp study was performed at 8 °C, which would significantly reduce ATP demand. To our knowledge, cardiac performance of H10 turtles during 2 h of anoxia is unsurpassed among vertebrates.

Similar to other studies [47, 52, 53], anoxic N21 and H10 snapping turtle hearts recovered full functionality when normoxia was restored. This rare adaptation is noteworthy, because atrial reoxygenation is nearly instantaneous in turtles [54], which could theoretically expose the heart to oxidative damage through the overproduction of reactive oxygen species (ROS). Although ROS are key signalling molecules and an inevitable product of electron transport, they can also be harmful by causing lipid peroxidation, protein and DNA damage, and apoptosis [55]. In mammals, ROS production during reperfusion is a major cause of ischemia/reperfusion injury [56]. In contrast, turtles fully recover contractile function after several hours of anoxia without any conspicuous injury [47, 57, 58]. The lack of reoxygenation injury in turtle hearts has been attributed to the maintenance of an ATP/ADP pool and low succinate accumulation, which reduces the likelihood of superoxide production [59]. Interestingly, we recently found that ROS production after anoxia was significantly lower in H10 vs. N21 snapping turtle cardiomyocytes [20], which may help to explain why H10 turtles in the present study could maintain higher levels of cardiac performance during reoxygenation (Fig. 1).

Effects of developmental hypoxia on cardiovascular regulatory mechanisms

Drawing from previous literature, we can make some inferences about the physiological and molecular mechanisms underlying differences in cardiac performance between H10 and N21 turtles. The H10 group had a blunted anoxic bradycardia compared to the N21 cohort, suggesting differences exist in factors that regulate heart rate, including local control mechanisms and/or intrinsic pacemaker properties. Acute anoxic bradycardia in warm turtles is mostly vagally mediated, but intrinsic pacemaker rate is also reduced by chronotropic extracellular factors such as acidosis, hyperkalemia, hypercalcemia, or adrenaline levels [46, 60]. Thus, developmental hypoxia could reduce the anoxia sensitivity of any of these pathways, leading to blunted anoxic bradycardia. Interestingly, H10 embryos are tachycardic when measured in normoxia and compared with N21 embryos. This is due to a blunted cholinergic tone, but the heart rate response to hypoxia was similar between groups [61]. We recently showed intracellular pH in anoxic snapping turtle cardiomyocytes is similar between N21 and H10 cohorts, but it is possible that extracellular pH, K⁺, Ca²⁺, or adrenaline levels differ, which could affect intrinsic pacemaker rate [60].

Enrichment analysis of differentially methylated genes in the current study revealed groups of genes that are functionally related to intrinsic pacemaker mechanisms, including GO terms “potassium channel complex,” “cation channel complex,” “transmembrane transporter complex,” and “ion channel complex.” We confirmed differential expression of one candidate from this category: sodium voltage-gated channel alpha subunit 5 (SCN5A) was up-regulated in hearts of 9-month-old turtles exposed to hypoxia as embryos. SCN5A mediates voltage-dependent sodium ion permeability of excitable membranes and inactivation of this channel is regulated by intracellular calcium. Knockout of SCN5A in mice results in embryonic lethality, with major defects in ventricular morphology [62]. More significantly, heterozygotes with one good copy of SCN5A survive, but display defects in atrioventricular and intramyocardial conduction, increased ventricular refractoriness, and tachycardia. Thus, elevated SCN5A expression could have the opposite effect and contribute to the higher heart rate of turtles from hypoxic incubations when exposed to acute anoxia.

Likewise, differentially methylated genes were enriched for genes that could alter autonomic control of heart rate, including GO Biological Process terms for “regulation of trans-synaptic signaling,” “regulation of nervous system development,” and “regulation of neuron differentiation,” and GO Cellular Component terms such as “synapse,” “glutamatergic synapse,” and “postsynaptic membrane.” In line with these observations, analysis of differentially methylated CGIs (i.e., potential enhancers) revealed significant enrichment of binding sites for Paired-Like Homeobox 2B (PHOX2B). This gene is intriguing, because it plays a key role in the formation of autonomic reflex pathways, including those controlling baroreflexes [63]. Mutation of PHOX2B in humans causes congenital central hypoventilation syndrome, which has cardiac arrhythmia as one of its incompletely penetrant symptoms [64]. Previous work has shown that hypoxic incubation alters α-adrenergic regulation of heart rate, β-adrenergic regulation of blood pressure, and influences expression of α- and β-adrenoreceptors in snapping turtle embryos [12]. Those findings point to hypoxic programming of the autonomic nervous system and/or tissue responsiveness to sympathoadrenal regulation. Further work is needed to identify mechanisms underlying the contribution of parasympathetic and sympathetic branches to the blunted bradycardic response in H10 turtles.

In addition to differences in the heart rate response to anoxia between experimental groups, anoxia led to an increase in cardiac inotropy in the H10 cohort, while this parameter was decreased in the N21 cohort. These results reflect the finding that N21 cardiomyocyte contractility remains depressed throughout 20 min of anoxia, while H10 contractility rebounds to pre-anoxic levels [20]. Improved anoxia tolerance of H10 cardiomyocytes was supported by enhanced myofilament calcium sensitivity and a superior ability to suppress ROS production [20]. Therefore, the increase in stroke volume observed in anoxic H10 turtles might be partly driven by intrinsic regulation of calcium and ROS. In support of this idea, we confirmed differential expression of peptidylprolyl isomerase A (PPIA) between N21 and H10 turtles. Secretion of this protein occurs in response to oxidative stress and plays a role in mediating Angiotensin II effects on cardiomyocyte hypertrophy by potentiating ROS production [65].

Enrichment analysis of differentially expressed genes in the current study identified other genes that might contribute more directly to positive inotropy in vivo and enhanced cardiomyocyte contractility in vitro: relevant GO Biological Process terms included “sarcomere organization” and “signal transduction,” as well as GO Cellular Component terms “z-disc” and “collagen trimer.” We confirmed differential expression of two candidate genes in this category: troponin T2, cardiac type (TNNT2) and tropomyosin 3 (TPM3) were both elevated in hearts of turtles exposed to hypoxia as embryos. Together, these proteins regulate calcium-dependent contraction of myofilaments. Increased TNNT2 and TPM3 expression could, therefore, cause differences in calcium sensitivity of cardiomyocytes in vitro [20] and contribute to the physiological differences observed in vivo in the present study.

Developmental hypoxia programs genome-wide DNA methylation and gene expression patterns

While examples of individual candidate genes discussed above are interesting, our study provides the first ever genome-wide analysis of DNA methylation and gene expression patterns in a reptile. WGBS allowed us to characterize the methylation landscape across the snapping turtle genome in unprecedented detail. Promoters and first exons displayed a bimodal pattern of CpG methylation, while other genome features displayed a unimodal distribution. CpGs were more highly methylated in genes than intergenic regions, which displayed a broader range of methylation levels. We also observed genome-wide correlations between CpG methylation and gene expression in the snapping turtle that are consistent with correlations observed in mammals, anuran amphibians, and fish [66]. In particular, methylation levels in promoters and first exons were negatively correlated with mRNA expression, whereas methylation in the remaining exons and introns was positively correlated with mRNA expression. Finer scale spatial analysis of CpG methylation across proximal promoters and the 5′ end of genes revealed a clear signature: higher methylation at TSSs was associated with lower expression, while lower methylation at TSSs was associated with higher expression (Fig. 5). These findings indicate DNA methylation near promoters plays a conserved role in repression of gene expression in turtles. The discovery of a positive correlation between methylation in gene bodies and gene expression in turtles is also observed in other vertebrate lineages [66]. It has been suggested that this positive correlation is a secondary effect of the greater accessibility of more highly transcribed genes to DNA methylating enzymes [67]. Together, these observations are significant, because DNA methylation is absent (or minimal) and plays no role in regulating gene expression in some model organisms [68, 69].

We also found that exposure to hypoxic conditions during embryogenesis programmed DNA methylation patterns and that methylation of CpGs vs. CGIs varied among genomic features. Intergenic regions, where enhancers and silencers are located, were enriched for differentially methylated CGIs but not for individual CpGs. If orphan CGIs in turtles act as enhancers, as suggested by recent studies in mammals [34,35,36], differential methylation of these sites might play an outsized role in driving differences in gene expression patterns between N21 and H10 turtles. In contrast to intergenic regions, promoters were less likely to contain differentially methylated CGIs but were enriched for differentially methylated CpGs. These findings are consistent with the observation that CGIs in promoters are usually unmethylated in mammals [37]. Finally, we found that exons displayed less differential methylation than expected by chance for both CpGs and CGIs. Overall, we detected significant relationships between hypoxia-induced DNA methylation, gene expression patterns, and cardiovascular physiology later in life, though links for individual gene are not linear and will require substantial experimental work to elucidate.

Nonetheless, we were able to gain insight into potential regulatory mechanisms through the identification of enriched sequence motifs (i.e., putative TF-binding sites) in promoters of differentially expressed genes. For instance, glucocorticoid response elements were enriched in proximal promoters of genes that were differentially expressed between ventricles of snapping turtles from hypoxic vs. normoxic incubations. Work in mammals shows that glucocorticoids play a role in maturation of the cardiovascular system during late gestation, including effects on peripheral resistance, blood pressure, and heart rate that protect against acute hypoxia in embryos [70]. Reciprocal interactions between hypoxia and glucocorticoids have also been observed: fetal hypoxia programs differential methylation and expression of the glucocorticoid receptor gene in rat hearts [71]. Although we did not detect differential methylation or expression of the glucocorticoid receptor gene in the snapping turtle, a potential link between hypoxia-induced differential gene expression and signaling via glucocorticoid response elements deserves further study. Perhaps there are programmed differences in the function of the hypothalamic–pituitary–adrenal axis between N21 and H10 turtles, as observed in adult rats exposed to intermittent hypoxia during the postnatal period [72].

Another proximal cis-regulatory element might play a major role in driving differential gene expression between hearts of N21 and H10 turtles. Promoters were enriched for sequence motifs recognized by Ras Responsive Element-Binding Protein 1 (RREB1). Mutations in RREB1 in humans cause Noonan-spectrum disorders, which are recapitulated in Rreb1 hemizygous mice [73]. Mice with one functional copy of Rreb1 exhibit cardiac hypertrophy and sensitization of cardiomyocytes to MAPK signaling. Interestingly, “signal transduction” was an over-represented GO term among differentially expressed genes, which included MAP3K5 and MAPK9, in turtle hearts. Anastasiadi et al. [66] also found enrichment of RREB1-binding sites in genes (gene body ± 4 kb) that display tissue-specific differential methylation. The observation that RREB1 is a methyl-CpG-binding TF suggests it may play a broader role in regulating differential gene expression in a DNA methylation-dependent manner [29].

We also identified sequence motifs that were significantly enriched within differentially methylated CGIs. We hypothesize these putative TF-binding sites are involved in the initial programming of differential DNA methylation of specific CGIs during embryogenesis as well as causing differential expression of target genes and differences in cardiovascular physiology later in life. Future studies using promoter capture HiC or HiCap could be used to test whether differentially methylated CGIs physically interact with promoters of differentially expressed genes [74, 75], as predicted by our model. However, it is not a straightforward task to link putative enhancers to their target promoters, because enhancers can act at distances of tens to hundreds of thousands of bases via DNA looping and skip over other promoters [74,75,76]. Those studies report that only a portion of enhancers (47% to 65%) physically interact with the nearest active TSS. Thus, the closest gene may not always be the target for putative “CGI enhancers” in the snapping turtle. This could explain why genes closest to differentially methylated CGIs were less likely rather than more likely to be differentially expressed in the snapping turtle. Another potential explanation is that the current version of the snapping genome is a scaffold level assembly. A fragmented genome could hamper our ability to physically associate differentially methylated CGIs with their closest targets.

Examination of transcription factors that bind enriched motifs in differentially methylated CGIs lends credence to our hypothesis. All three hypoxia inducible transcription factors (HIF1A, ARNT, and EPAS1) are in the list of enriched motifs and are candidates for programming differential DNA methylation in turtle hearts (Additional file 4: Table S4). Hypoxia can cause global changes in DNA methylation by regulating expression of DNA methyltransferases and ten–eleven translocation (TET) methylcytosine dioxygenases [77,78,79]. Yet, the hypoxia-induced methylation patterns observed here are specific rather than global and suggest targeting to cis-regulatory elements of genes that are differentially expressed between turtles from hypoxic vs. normoxic incubations. In fact, recent work shows that DNA methylation directly interferes with HIF binding and that cell-type specific methylation patterns determine responsiveness of HIF target genes to acute hypoxia [79]. We suggest an analogous mechanism could underlie programmed differences in DNA methylation and gene expression patterns between N21 and H10 turtles. Activation and binding of HIFs and TFs such as RREB1 and SMADs to specific sites would drive differential methylation of CpGs and CGIs in embryonic hearts. These programmed patterns would then cause differential gene expression and differences in cardiovascular phenotype later in life. Differentially methylated CGIs in turtle ventricles were enriched in binding sites for transcription factors associated with cardiovascular, mitochondrial, or autonomic defects in mammals (e.g., CUX2, GABPA, GSC, PHOX2B, SMAD4, and ZEB2).

SMAD4-binding sites were enriched within differentially methylated CGIs in turtle ventricles. Given that SMAD4 binds methylated sites [29], the intersection of DNA methylation and TGF-β signaling in the turtle heart is particularly intriguing. Indeed, embryonic hypoxia programmed differential expression of BMP10 and BMPR2 in turtle hearts. BMP10 is a key signaling molecule in the developing and mature heart in mammals [80, 81]. BMP10 is a member of the TGF-β protein family and binds to BMPR2 and ALK1 to trigger phosphorylation of SMAD2/3 or SMAD1/5/8, which form complexes with SMAD4. SMAD complexes then translocate to the nucleus to regulate expression of target genes. Signaling via SMAD4 in cardiomyocytes plays a crucial role in regulating sarcomere function, ion-channel gene expression, cardiomyocyte survival, and cardiac function in adult mice [82]. Gain of function mutations in SMAD4 in humans cause Myhre Syndrome with cardiomyopathy [83]. Mutations of BMPR2 in humans and mice also cause cardiovascular phenotypes that are associated with pulmonary arterial hypertension [84, 85]. Hautefort et al. [86] have recently shown that rats heterozygous for a BMPR2 exon deletion display changes in right ventricular cardiomyocyte morphology and physiology, including smaller diameter, decreased calcium sensitivity, and decreased contractility. Thus, hypoxia-induced changes in TGF-β signaling via the BMP10/BMPR2/SMAD4 axis could direct differential DNA methylation during embryogenesis and influence subsequent cardiac physiology in N21 and H10 turtles. In fact, TGF-β-dependent activation and binding of SMADs has been shown to displace the ZNF217/CoREST/DNMT3A complex and cause demethylation of the p15^ink4b promoter and induce expression of the p15^ink4b gene [87]. There is also potential for crosstalk between TGF-β signaling and MAP kinase signaling, because SMADs and RREB1 (discussed above) have been shown to directly interact to regulate epithelial to mesenchymal transitions and induce fibrosis in myofibroblasts [88].

Conclusion

Overall, our study shows chronic hypoxia during embryogenesis significantly improves cardiac anoxia tolerance in juvenile snapping turtles and that these effects are associated with changes in DNA methylation and gene expression patterns. Our findings also point to specific genes and signaling pathways that may underlie extreme hypoxia/anoxia tolerance in the snapping turtle. Based on these findings, we propose a model in which hypoxia during embryogenesis activates hypoxia inducible factors (HIF1A, ARNT, and EPAS1) and other key TFs (e.g., RREB1 and SMAD4), which interact with specific-binding sites to direct (or inhibit) methylation of nearby CpGs and CGIs (Fig. 8). Hypoxia-induced DNA methylation patterns would then be passed down through cell divisions and maintained in later life (i.e., they are programmed). Differential methylation of CpGs and CGIs would modulate promoter/enhancer/silencer activity, chromatin structure, and influence gene expression by affecting binding of the same transcription factors (e.g., HIFs, RREB1, and SMAD4) later in life. This, in turn, would drive differences in cardiomyocyte and cardiac physiology. It is important to point out that these findings are correlative and that further research will be required to test the hypothesized mechanisms. There are limitations in the approaches that can be used in turtles. Genetic approaches such as gene knockouts and transgenics are not yet feasible with turtles. However, potential experimental approaches include primary cell culture of cardiomyocytes. An in vitro system would be amenable to pharmacological manipulations of the signaling pathways identified here. Mechanisms could also be studied by transient transfection of expression vectors with candidate genes from this study and/or knockdown of genes with siRNA or lentiviral vectors carrying shRNA.

Fig. 8

Hypothetical model for developmental programming of DNA methylation and gene patterns by hypoxic incubation. Exposure to low oxygen during embryogenesis activates HIFs and other transcription factors (TFs), which bind to specific sites in the genome of developing cardiomyocytes. These factors recruit DNMTs and TET in a locus specific manner to methylate and demethylate adjacent CpGs, respectively. DNA methylation patterns in CpGs and CGIs are inherited mitotically and influence gene expression patterns later in life

The physiological significance of these findings also awaits further research. On one hand, increasing cardiac output during acute anoxia might be beneficial for breath-hold dives when snapping turtles are foraging. During overwintering periods, this increased capacity, combined with the low ambient water temperature, could allow the animals to sustain activity if needed. However, this strategy could become risky for longer periods of anoxia if ATP turnover is elevated and glycogen reserves become limited. In this regard, it would be interesting to measure levels of lactate production in anoxic H10 turtles to assess ATP turnover rate.

Methods

Turtle collection, incubation, and husbandry

Common snapping turtle (Chelydra serpentina) eggs were collected from the wild in Minnesota, USA, and transported to the University of North Texas, TX, USA, for incubation. Permission to collect the eggs was granted to DA Crossley by the Minnesota Department of Natural Resources (permit no. 21232). On arrival two eggs from individual clutches were staged to establish embryonic age. Eggs were embedded to their midpoint in vermiculite, inside plastic boxes (2.5-L Ziploc® containers, SC Johnson, Racine, WI, USA) that were placed inside large (75.7 L) sealable plastic bags (Ziploc®). All incubation conditions were carried out in a walk-in environmental control room (model IR-912L5; Percival Scientific, Perry, IA, USA). The vermiculate was mixed in a 1:1 ratio with water, as previously described [89]. Incubation lasted no more than 55 days and all eggs were maintained at 30 °C, a female-determining temperature [61].

At approximately 20% of development (9–12 days after laying; determined by embryonic staging), eggs were randomly assigned to either normoxic/atmospheric oxygen (21% O₂; designated as N21) or hypoxic (10% O₂; designated as H10) cohorts for the remainder of embryonic development. To achieve the desired oxygen level, parallel gas inflow and outflow tubes were attached to the large Ziploc® bags and O₂ gas mixtures were set at a flow rate of 2–3 L min⁻¹ using rotameters (Sho-Rate Brooks Instruments Division, Hatfield, PA, USA) downstream of either compressed N₂/air mixture or air alone. The gas mixtures passed through an H₂O bubbler to ensure 80–95% relative humidity and their compositions were monitored continuously with an oxygen analyser (S-3AI; AEI Technologies, Pittsburgh, PA, USA).

Upon hatching, all turtles were housed in common, normoxic (21% O₂) conditions at 26 °C, in separate normoxic and hypoxic groups. Turtles were fed ad libitum, with dry turtle food (Mazuri, PMI Nutrition International, Brentwood, MO, USA) 2–4 times weekly and kept in a daily 12:12 light–dark cycle until experimentation (up to 1.5 years).

In situ turtle cardiac anoxia tolerance

Size- and clutch-matched turtles from each developmental cohort were studied 1.5 years after hatching (N = 6 and 5, for N21 and H10, respectively). Turtle body and heart masses are provided in Table 1. Prior to experimentation, turtles were anaesthetized in a sealed box containing cotton gauze saturated in isoflurane (Isoflo®, Abbott Laboratories, North Chicago, IL, USA). Once pedal and eye reflexes were absent, turtles were removed from the box, placed ventral-side up and intubated with flexible Tygon® tubing that was inserted into the trachea via the glottis. A ventilator (model 683, Harvard Apparatus, Holliston, MA, USA) and vaporizer (FluTec vaporizer, FluTec, Ohmeda, OH, USA) provided mechanical ventilation with 3% isofluorane, at a rate of 3–4 breaths min⁻¹ and tidal volume of 20 mL kg⁻¹. A gas-mixer (GF-3mp, Cameron Instrument Company, Port Aransas, TX, USA) was connected to the ventilator and controlled the composition of gases.

A square cut (4 cm²) was made in the plastron directly over the heart to expose the major cardiac outflow vessels and pericardium. Major arteries were isolated from surrounding tissue by blunt dissection for placement of the blood-flow probes (Transonic Systems, Ithica, NY, USA). One probe (3- or 4-mm diameter) was used to measure blood flow in the right aorta, both subclavian arteries, and the right carotid collectively. Separate probes were used to measure blood flow in the left aortic, left carotid artery (both 1–2 mm), and the left pulmonary artery (1.5–2.5 mm). Each flow probe was calibrated at 30 °C, with an infusion syringe pump (PHD 2000, Harvard Apparatus, USA). The flow probes were connected to two T206 blood-flow meters (Transonic Systems Ithica, NY, USA). To measure ventricular pressure, a small hole was made in the apex of the heart using a 22-gauge needle and a pressure catheter (size 1.4 F, model SPR-671, Millar Instruments, Houston, TX, USA) was inserted into the lumen of the heart. The catheter was connected to an amplifier (MPVS-300, Millar Instruments) which was calibrated daily against a static column of water, using a two-point calibration (0 and 1 kPa). The outputs from the flowmeters and pressure amplifier were connected to a PowerLab® 8/35 data-recording system (ADInstruments, Colorado Springs, CO, USA) and recorded on a computer, with LabChart Pro® software (v8.2, ADInstruments), and data were recorded at 100 Hz.

After the flow probes and catheters were placed, isofluorane was reduced to 1–1.5%, ventilation was raised to 10–11 breaths min⁻¹, and cardiovascular variables were left to stabilize for at least 30 min before the experimental protocol commenced. The experiment was designed to measure cardiac function during three distinct periods: 10 min of normoxia (21% O₂, 3% CO₂, and 76% N₂), 120 min of anoxia (3% CO₂ and 97% N₂), and 30 min of reoxygenation (21% O₂, 3% CO₂, and 76% N₂). The ventilated gas mixture was regularly checked with oxygen and carbon-dioxide analyzers (model S-3A/I and CD-3A, respectively, Ametek, Berwyn, PA, USA). All studies were carried out according to an approved animal-care protocol of the University of North Texas Institutional Animal Care and Use Committee (no. 1403-04).

Mean blood-flow (Q̇) values were calculated from the average of 5-min data periods throughout the experimental protocol. Total systemic blood flow (Q̇_Sys) was calculated as the sum of flow from the right and left aortas, subclavian arteries, and carotid arteries, whereas total pulmonary blood flow (Q̇_Pul) was calculated as 2× the flow of the left pulmonary artery, assuming that flows through the left and right pulmonary arteries are identical. Total cardiac output (Q̇_Tot) was calculated as the sum of Q̇_Sys and Q̇_Pul. Total, systemic, and pulmonary stroke volumes ($V_{{{\text{S}},{\text{Tot}}}}$, $V_{{{\text{S}},{\text{Sys}}}}$, and $V_{{{\text{S}},{\text{Pul}}}}$, respectively) were calculated using the following equation:

$$ V_{{\text{S}}} = \frac{{\dot{Q}}}{{f_{{\text{H}}} }}, $$

(1)

where Q̇_Tot, Q̇_Sys, and Q̇_Pul, were used to find $V_{{{\text{S}},{\text{Tot}}}}$, $V_{{{\text{S}},{\text{Sys}}}}$, and $V_{{{\text{S}},{\text{Pul}}}}$, respectively.

Net and fractional shunts were calculated using Eqs. 2 and 3, respectively, to assess the distribution of blood flow between the pulmonary and systemic circulations:

$$ \dot{Q}_{{{\text{Shunt}}}} = \dot{Q}_{{{\text{Pul}}}} - \dot{Q}_{{{\text{Sys}}}} , $$

(2)

$$ \dot{Q}_{{{\text{Fractional}}}} = \frac{{\dot{Q}_{{{\text{Pul}}}} }}{{\dot{Q}_{{{\text{Sys}}}} }}. $$

(3)

Mean ventricular pressure ($P_{{{\text{Vent}}}}$) was calculated using the following equation:

$$ P_{{{\text{Vent}}}} = \frac{{P_{{{\text{Systolic}} }} + 2P_{{{\text{Diastolic}}}} }}{3}, $$

(4)

where $P_{{{\text{Systolic}} }} \;{\text{and}}\;P_{{{\text{Diastolic}}}}$ are systolic and diastolic pressure, respectively.

Finally, cardiac power output (PO) was calculated using the following equation:

$$ {\text{PO }} = \frac{{\dot{Q}_{{{\text{Total}}}} \cdot \Delta P}}{{{\text{heart}}\;{\text{mass}}}}, $$

(5)

where $\Delta P = P_{{{\text{Systolic}} }} - P_{{{\text{Diastolic}}}}$.

Transcriptome analysis of hearts exposed to developmental hypoxia

RNA-Sequencing (RNA-Seq) was carried out to measure steady state differences in cardiac gene expression between juvenile turtles exposed to normoxic or hypoxic conditions during embryogenesis. Normoxic and hypoxic groups included 7-month-old (n = 5) and 9-month-old turtles (n = 3), for a fully factorial design (total n = 16). The hypoxic group included equal numbers of turtles with normal-sized (n = 4) and enlarged hearts (n = 4) relative to their body size.

Hearts were dissected from turtles and weighed. Atria and ventricles were separated, placed in microfuge tubes, snap frozen in liquid nitrogen, and stored at − 80 °C. Total RNA was isolated from ventricles by grinding frozen tissue with a mortar and pestle on dry ice. Frozen, pulverized tissue was transferred to a tube containing Trizol and homogenized for another 30 s using a BioGen PRO200 homogenizer with a 5 mm generator probe. Remaining steps were carried out according to the manufacturer’s protocol. The only modification was 2 additional extractions with 500 µL of chloroform to remove phenol traces from the aqueous phase prior to RNA precipitation. RNA quality was high with RINs ranging from 8.4 to 9.1 and no indication of genomic DNA contamination when assessed on agarose gels or via qPCR (no amplification).

Total RNA was used as input for the NEB PolyA nondirectional library preparation kit. Barcoded cDNA libraries with 250–300 bp insert sizes were sequenced on Illumina HiSeq system (150 bp, paired end reads) by Novogene. One set of 6 samples was sequenced to a depth of 29 to 33 million raw reads (forward + reverse), while a second set of 10 samples was sequenced to a depth of 65 to 150 million raw reads (forward + reverse) (Additional file 5: Table S5). The first 11 bp of reads were cropped and low-quality bases trimmed (sliding window of 4 bp and average Q score ≥ 15) with a minimum read length of 30 bp using Trimmomatic [90]. Reads were mapped to the snapping turtle genome [91] using HISAT2 with default parameters [92]. featureCounts [93] was used to extract read counts from BAM files for subsequent gene expression analyses.

DESeq2 was used to screen for differences in gene expression [94]. Oxygen concentration, age, and the oxygen concentration by age interaction were independent factors in a general linear model with a genewise P < 0.01. The distribution of FPKMs were manually examined to identify differences driven by outliers. Two-way ANOVAs were then carried out on FPKMs for each gene identified by DESeq2 to ensure differences were significant using two different statistical models. DESeq2 uses a hierarchical model with likelihood ratio tests and shrinks estimates of dispersion by assuming genes with similar expression values display similar variance. In contrast, ANOVA employs ordinary least squares with F tests that use empirically derived variance estimates for each gene. Genes were excluded from the final list of differentially expressed genes when an outlier drove a significant effect or when DESeq2 and two-way ANOVA results were not concordant (i.e., differences were not robust to the statistical model). Gene expression was also compared between normal-sized hearts (n = 12) vs. enlarged hearts (n = 4) with an FDR adjusted p value < 0.1. The final set of differentially expressed genes included those affected by oxygen concentration, the oxygen concentration by age interaction, and the genes that differed between normal-sized (both N21 and H10) vs. enlarged hearts (H10). Genes that only changed with age were not analyzed any further, because the long-term effect of hypoxia was the primary focus of this study.

Validation of differential gene expression in hearts exposed to developmental hypoxia

qPCR was used to measure gene expression in a larger set of samples from the same experiment that produced animals for the RNA-Seq study (i.e., 13 normoxic hearts and 12 hypoxic hearts; 20 normal-sized hearts and 5 enlarged hearts). Total RNA was extracted as described above. Reverse transcription and absolute qPCR with rigorous standard curves were carried out as previously described [61, 95]. Expression of CACNA2D1, CNP, and YTHDF3 were not affected by any independent variables so these genes were used as controls. The first component from a principal components analysis of these genes was used as a covariate for analysis of the remaining genes. This covariate serves as a control for variation in the quality of input RNA and the efficiency of reverse transcription reactions (i.e., such as a housekeeping gene).

Methylome analysis of hearts exposed to developmental hypoxia

Turtles exposed to normoxic (n = 3) or hypoxic conditions (n = 3) during embryogenesis were used for WGBS. DNA was extracted from frozen, pulverized ventricles of the same 9-month-old turtles used for the RNA-Seq study. DNA was extracted using the DNeasy Blood and Tissue kit from Qiagen. Agarose gel electrophoresis of DNA revealed high molecular weight DNA (> 60 kb) with no RNA contamination. Six μg of DNA was shipped to Novogene for WGBS. Libraries were prepared with 200–400 bp insert sizes. Bisulfite (BS) conversion was carried out with the EZ DNA Methylation Gold Kit from ZymoResearch. Libraries were sequenced on a NovaSeq 6000 instrument. QC analysis of raw reads showed BS conversion rate was greater than 99.9% for all libraries and coverage ranged from 32.9× to 37.6× (Additional file 6: Table S6).

Trimmomatic was used to remove TruSeq3 PE adapters, trim 3 bp from the 5′ and 3′ ends, and trim low-quality bases (sliding window of 4 bp and average Q score ≥ 15) with a minimum read length of 36 bp [90]. Reads were mapped to the snapping turtle genome using the Bismark bisulfite read mapper [96]. Mapping statistics for each library are summarized in Additional file 7: Table S7. Average mapping efficiency was 79.2%, which is excellent for WGBS data [97]. As expected, methylated cytosines were primarily found in the context of CpG dinucleotides (75%). Few methylated cytosines were found in the context of CHG (0.2%) or CHH (0.2%) trinucleotides, where H is any base except G. Approximately 3.4% of methylated cytosines were in an unknown context.

methylKit was used to call methylated CpGs and determine whether methylation levels were significantly different between N21 and H10 groups [98]. A minimum coverage of 10 in two of three replicates was required for statistical comparison. Differences between N21 and H10 groups were called significant for individual CpGs if the difference in methylation was > 25% and q < 0.01 (q is the FDR adjusted p value).

We used the newcpgreport tool (https://www.bioinformatics.nl/cgi-bin/emboss/newcpgreport) to call CGIs in the snapping turtle genome using default parameters: Obs/Exp > 0.6, %C + %G > 50, and length > 200 bp. We identified 201,828 CGIs in the snapping turtle genome which is less than the 307,193 CGIs in the human genome with the same parameters [37]. When corrected for genome size, however, the frequency of CGIs is similar at 89,383 CGIs/Gb in the snapping turtle and 93,089 CGIs/Gb in humans. Overall methylation of CGIs was calculated as the sum of methylated CpGs divided by the total number of CpGs within an island, which is essentially the average % methylation across the island. Comparisons between N21 and H10 groups were made using the Fisher Exact test and q < 0.05.

Statistical analyses

Data were analyzed for statistical significance by a mixed-effects, generalized linear model (GLM), using Šidák post-hoc corrections for pairwise comparisons, with SPSS 25 (IBM, Armonk, NY, USA). For the GLMs, developmental oxygen (normoxia or hypoxia), acute oxygen treatment (normoxia, anoxia, or reoxygenation), and time were the independent variables and cardiovascular variables were the dependent variables. Significance was accepted when p ≤ 0.05. All data are reported as means ± standard error (SEM). GOATOOLS [99] was used to test for functional enrichment of gene ontology (GO) terms among differentially expressed genes from the RNA-Seq study and differentially methylated genes from the WGBS study. Genes identified in those experiments were compared to a species-specific list of GO terms generated by Das et al. [91].