Transcriptomic analysis of preovipositional embryonic arrest in a nonsquamate reptile (Chelonia mydas).

After gastrulation, oviductal hypoxia maintains turtle embryos in an arrested state prior to oviposition. Subsequent exposure to atmospheric oxygen upon oviposition initiates recommencement of embryonic development. Arrest can be artificially extended for several days after oviposition by incubation of the egg under hypoxic conditions, with development recommencing in an apparently normal fashion after subsequent exposure to normoxia. To examine the transcriptomic events associated with embryonic arrest in green sea turtles (Chelonia mydas), RNA-sequencing analysis was performed on embryos from freshly laid eggs and eggs incubated in either normoxia (oxygen tension ~159 mmHg) or hypoxia (<8 mmHg) for 36 h after oviposition (n = 5 per group). The patterns of gene expression differed markedly among the three experimental groups. Normal embryonic development in normoxia was associated with upregulation of genes involved in DNA replication, the cell cycle, and mitosis, but these genes were commonly downregulated after incubation in hypoxia. Many target genes of hypoxia inducible factors, including the gene encoding insulin-like growth factor binding protein 1 (igfbp1), were downregulated by normoxic incubation but upregulated by incubation in hypoxia. Notably, some of the transcriptomic effects of hypoxia in green turtle embryos resembled those reported to be associated with hypoxia-induced embryonic arrest in diverse taxa, including the nematode Caenorhabditis elegans and zebrafish (Danio rerio). Hypoxia-induced preovipositional embryonic arrest appears to be a unique adaptation of turtles. However, our findings accord with the proposition that the mechanisms underlying hypoxia-induced embryonic arrest per se are highly conserved across diverse taxa.

INTRODUCTION

Hypoxia can induce suspended animation in diverse eukaryotic taxa, including yeast (e.g., Saccharomyces cerevisiae Chan & Roth, 2008) and the embryos of some nematodes (e.g., Caenorhabditis elegans, Miller & Roth, 2009), insects, (e.g., Drosophila melanogaster, DiGregorio et al., 2001), and fish (e.g., Danio rerio, Padilla & Roth, 2001). A number of genes and associated mechanisms have been identified as important mediators of developmental arrest in these taxa. For example, in the yeast Saccharomyces cerevisiae, genes with roles in ribosome biogenesis, transcription, translation, and many associated processes were significantly downregulated during arrest induced by hypoxia, as were genes involved in mitochondrial retrograde signalling (Chan & Roth, 2008). In embryos of the nematode Caenorhabditis elegans, san1, encoding senataxin‐associated nuclease 1 (SAN‐1), a component of the spindle assembly checkpoint, was found to inhibit progression from the M to G phase of the cell cycle during hypoxia (Miller & Roth, 2009; Nystul et al., 2003). In embryonic zebrafish (Danio rerio), hypoxia‐induced arrest appears to be mediated by hypoxia‐induced upregulation of insulin‐like growth factor binding protein‐1 (IGFBP‐1) via a hypoxia inducible factor‐1 (HIF‐1) dependent mechanism, resulting in inhibition of IGF‐1 and IGF‐2 stimulated cell proliferation (Kajimura et al., 2005, 2006). Nitric oxide appears to be a critical factor in suspended animation and hypoxia tolerance in Drosophila embryos (Teodoro & O'Farrell, 2003).

Embryos of both freshwater and marine turtles enter a period of arrested development after gastrulation and morphological development of the embryo does not recommence until after oviposition (Rafferty & Reina, 2012). This preovipositional embryonic arrest appears to be induced and maintained by the hypoxic environment of the turtle oviduct (Kennett et al., 1993; Rafferty et al., 2013; Rings et al., 2015; Williamson et al., 2019; Williamson, Evans, & Reina, 2017; Williamson, Evans, Robinson, & Reina, 2017). After eggs are laid, exposure to atmospheric levels of oxygen triggers recommencement of development (Williamson, Evans, & Reina, 2017). Preovipositional embryonic arrest appears to be an adaptation that enables female turtles to have greater flexibility in their reproductive schedule. Further, it provides a survival benefit for offspring due to synchronization of post‐ovipositional development and thus their time of hatching, as a predator‐swamping strategy (Santos et al., 2016). Preovipositional arrest is not apparent in other oviparous reptiles such as crocodilians (Williamson, Evans, Manolis, et al., 2017), but has been observed in all turtle species in which it has been studied to date (Kennett et al., 1993; Rafferty et al., 2013; Rings et al., 2015; Williamson et al., 2019; Williamson, Evans, & Reina, 2017; Williamson, Evans, Robinson, & Reina, 2017).

Other than an appreciation for the role of oxygen in triggering the breaking of preovipositional embryonic arrest in turtles, little is known of the mechanisms that underlie this form of embryonic arrest at the molecular and cellular level. Indeed, while it is well established that hypoxia arrests detectable morphological development, its effect on transcriptional activity remains to be determined. Embryonic arrest can be prolonged after oviposition by incubation of freshly‐laid turtle eggs in hypoxia (Kennett et al., 1993; Rafferty et al., 2013; Rings et al., 2015; Williamson et al., 2019; Williamson, Evans, & Reina, 2017; Williamson, Evans, Robinson, & Reina, 2017). In at least some species of turtle, prolongation of embryonic arrest after oviposition results in reduced hatching success, mainly due to early embryonic mortality (Rafferty et al., 2013; Rings et al., 2015) and may even affect morphology and fitness of hatchlings (Rings et al., 2015). These observations are consistent with the proposition that this form of embryonic arrest is a dynamic process and thus should be associated with observable changes in the transcriptome.

In the current study, we compared the transcriptome of embryos from freshly laid eggs of the green sea turtle (Chelonia mydas) with those of eggs that had been incubated for 36 h in either hypoxic (oxygen tension <8 mmHg) or normoxic (~159 mmHg) conditions. Thus, our experimental design provided an opportunity to gain information on the transcriptional events that normally occur with the breaking of embryonic arrest when turtle eggs are exposed to atmospheric levels of oxygen, and when arrest is extended by incubation in hypoxia. We observed profound differences in the transcriptome of embryos after extended arrest compared with embryos of freshly‐laid eggs, indicating the dynamic nature of this form of suspended animation. We also found patterns of gene expression in turtle embryos consistent with mechanisms of hypoxia‐induced suspended animation in other taxa, consistent with the proposition that these mechanisms are highly conserved.

METHODS

Permits and animal ethics

All experimental procedures were approved by the Animal Ethics Committee of Monash University's School of Biological Sciences (approval BSCI/2018/15), in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes. The research was conducted under a scientific permit issued by the Queensland Department of Environment and Heritage Protection (PTU18‐001205).

Egg collection and experimental design

Oxygen consumption of green sea turtle eggs is extremely low for the first 20–25 days of incubation, at least in comparison to avian eggs, reflecting the slow development of turtle embryos (Ackerman, 1980). A consequence of this low level of oxygen consumption is that the partial pressure of oxygen inside newly laid clutches of sea turtle eggs approaches that of atmospheric oxygen concentration at sea level (~159 mmHg; Ralph et al., 2005; Wallace et al., 2004). In contrast, the partial pressure of oxygen within the oviduct of gravid green sea turtles when they come onto the beach to oviposit is <5 mmHg (Rafferty et al., 2013). Our experimental treatments were designed to simulate the physiological conditions within the nest after oviposition (36 h incubation under normoxic conditions) and during extension of arrest within the oviduct (36 h incubation under hypoxic conditions) so embryos incubated under these conditions could be compared both with each other and with embryos from freshly laid eggs (0 h). The 36‐h time‐point was chosen because, in eggs incubated under normoxic conditions, it represents a stage of early embryonic development after embryonic arrest has been irreversibly broken, which in green sea turtle occurs 12–16 h after oviposition (Williamson, Evans, & Reina, 2017). Thus, comparison of the transcriptome of embryos incubated in normoxia for 36 h relative to that of embryos from freshly laid eggs has the potential to provide insight into the processes that normally occur after the breaking of preovipositional embryonic arrest in sea turtles. Conversely, comparison of the transcriptome of embryos incubated in hypoxia for 36 h with those from embryos in the other two treatments has the potential to provide insight into the mechanisms that maintain embryonic arrest, both before oviposition and when arrest is artificially maintained after oviposition.

Three eggs were collected from each of five female green turtles (Chelonia mydas) at Heron Island, Great Barrier Reef, Australia (23°26′18.71″ S, 151°54′30.23″ E) during the nights of 5–6 November 2018. Eggs were collected into plastic bags as they were laid and rapidly transported to the on‐site laboratory (within 15 min after collection). Breaking of arrest in green sea turtle embryos does not occur until 12–16 after oviposition (Williamson, Evans, & Reina, 2017). Thus, this short period of contact with atmospheric levels of oxygen during collection and transport of eggs is unlikely to have greatly confounded our findings. Three experimental conditions were then established: incubation for 36 h in normoxia (159 mmHg, 36 h N) before removal of the embryo (n = 5), incubation for 36 h in hypoxia (<8 mmHg, 36 h H) before removal of the embryo (n = 5), or immediate removal of the embryo (0 h, n = 5). One egg from each clutch (i.e., each female) was assigned to each group.

Embryos from the 0 h group were isolated immediately after arriving at the laboratory (see details below), within 1 h of oviposition. Eggs from the 36 h N group were buried in sand in a polystyrene container and incubated under normoxic conditions at 26–28°C in the laboratory. Eggs for the 36 h H group were incubated under hypoxic conditions in a sealed plastic bag filled with nitrogen (Pac Food Pty Ltd), as described previously (Williamson, Evans, Robinson, & Reina, 2017). This bag was in turn placed in a larger bag which was also filled with nitrogen gas and sealed. The bag was then buried in sand within a polystyrene container for incubation at 26–28°C in the laboratory. Embryos were removed after 36 h of incubation in normoxia or hypoxia.

Embryo isolation

The eggs were cleaned of adhering sand and debris and placed on a petri dish, stabilized with an aluminium foil ring. The method of isolation, which differs slightly for eggs before and after development of the white spot, has been described in detail elsewhere (Garriz et al., 2020). Once isolated, embryos were immediately frozen in liquid nitrogen before being stored at −80°C. Samples were then transported, on dry ice (−78.5°C), from Heron Island to Monash University by boat, airplane and car (total trip duration ~13 h) to a storage facility where they were again held at −80°C until further analysis.

Total RNA extraction and sequencing

RNA was isolated using QIAzol lysis buffer (Qiagen) following the manufacturer's protocol. The aqueous fraction was purified using the RNeasy Micro Isolation Kit (Qiagen), and RNA was eluted in 12 μl of RNase‐free distilled water. RNA samples were stored at −80°C until further analysis.

The concentration and purity of RNA obtained from embryos was assessed using capillary electrophoresis (CE) integrity analysis on an Advanced Analytical Technologies fragment analyser (Agilent Technologies) with fluorometric quantitation (Invitrogen Qubit). RNA quality was determined by examination of RNA size distribution on RNA Nano LabChips (Agilent) processed on a bioanalyser (Agilent 2100). RNA quantity was at least 1 ng/μl for each sample and the RNA integrity number (RIN; Schroeder et al., 2006) was ≥7 for all samples, so considered viable for RNA‐sequencing analyses (Jahn et al., 2008).

RNA samples were submitted to Micromon Genomics (Monash University) for library construction using an Illumina TruSeq RNA Sample Prep Kit version 2 (Illumina Inc.) following the manufacturer's instructions. Samples were subjected to 18 cycles for library amplification with 100 ng of RNA. The complete and prepared 76 bp SE libraries were sequenced using an Illumina Next‐Seq500 sequencing by synthesis (SBS) with SBS chemistry version 1.5 prepared according to the manufacturer's instructions (Illumina document #15046563 version 06; Illumina Inc.). Raw data were then generated with bcl2fastq2 conversion software (Illumina Inc.).

Data analysis and bioinformatics

The sequencing data were processed using the Laxy tool (Perry & Powell, 2020) and the RNA‐Sik version 1.5.4 pipeline (Tsyganov et al., 2018) using default parameters against the Chelonia mydas reference genome rCheMyd1 (GCF_015237465.1, NCBI). This included mapping the reads using the spliced transcripts alignment to a reference (STAR) RNA‐sequencing aligner version 2.7.2b (Dobin et al., 2013). Soft clipping of reads was performed by STAR as a default, with quantification against the release version 101 annotation using the featureCounts version 1.6.2 program (Liao et al., 2014). This resulted in an alignment rate to the reference of 83%–93% and gene assignment of approximately 65%, depending on the sample. Raw counts were analysed with Degust version 4.1 (Powell et al., 2019), a web tool which performs differential expression analysis using the Limma‐Voom package (limma_3.40.6) in the R version 3.6.1 statistical software environment (Law et al., 2014; Ritchie et al., 2015). Raw read counts were normalized by two methods: (i) counts per million (CPM) library size normalization and (ii) trimmed mean of M values (TMM) normalization (Robinson & Oshlack, 2010). Differentially expressed genes were defined as those with a false‐discovery rate (FDR) of ≤0.05, showing a >2‐fold (log2FC ≥ ±1) change in expression. A heatmap was generated by Heatmapper using unsupervised cluster differential expression analysis (Babicki et al., 2016). The RNA sequence data have been deposited in the NCBI Gene Expression Omnibus (Edgar et al., 2002) with accession no. GSE197628 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE197628).

Functional enrichment and pathway analysis

Over‐representation analysis (ORA) and gene set enrichment analysis (GSEA) methods were used to perform enrichment and pathway analysis. For ORA, lists for positive and negative fold changes were generated for all three possible between‐group comparisons (i.e., 0 h vs. 36 h N; 0 h vs. 36 h H, and 36 h H vs. 36 h N). Enrichment analysis with the ORA method was performed using Gprofiler (Raudvere et al., 2019). Few green sea turtle gene functions have been validated in the organism itself and no gene ontology (GO) annotations are available. Thus, orthologue analysis was performed by mapping genes with high sequence similarity and transferring functional annotations available for the freshwater painted turtle (Chrysemys picta bellii). GO, Reactome and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used for enrichment analysis. p‐values adjusted for false discovery rate (Benjamini‐Hochberg) with a cutoff of < .05 were used to filter altered GO terms and pathways in each pairwise comparison. Redundant parent GO terms were removed using OntologyIndex in the R package to focus on more fine‐level GO results. The results for the top five enriched terms for each comparison were combined and figures were generated using the R Statistical package (R Core Team, 2021). GSEA was carried out using fgsea in the R statistical software environment using all the observed genes in the RNA‐sequence data set. The genes in each pairwise comparison were ranked based on their fold‐change and false discovery rate‐adjusted P values. The gmt file containing multiple pathways, downloaded from Reactome.org, was used as a gene set list. The top up‐ and downregulated pathways were filtered based on normalized enrichment score (NES) and were visualized using ggplot in the R statistical software environment.

RESULTS

Morphological observations

The absence or presence of development of the embryo was assessed by visual inspection of the eggs. The characteristic white spot on the upper surface of the egg is indicative of the breaking of arrest, since it occurs once the vitelline membrane attaches to the shell membrane at the top of the egg, which in turn only occurs once the embryo develops past the gastrula stage and neurulation begins (Thompson, 1985). All five eggs incubated in normoxia (36 h N) developed a “white spot” (see Figure 1 in Garriz et al., 2020). In contrast, none of the five eggs incubated in hypoxia (36 h H) developed a white spot, indicating lack of development of the embryo past the gastrula stage. The absence of morphological development during hypoxic incubation, past Yntema's stage 0 (Yntema, 1968) and Miller's stage 6 (Miller, 1985) has previously been confirmed for Chelonia mydas (Rafferty et al., 2013; Williamson, Evans, & Reina, 2017) and other turtle species (Rafferty et al., 2013).

FIGURE 1

Heat map showing the relative expression of the 250 transcripts that varied most between the three conditions (p < .0001). Embryos were removed from eggs immediately after they were transported to the laboratory (0 h, embryos 1–5) or after incubation of the eggs at 26–28°C for 36 h in normoxia (36 h normoxia, embryos 6–10) or 36 h in hypoxia (36 h hypoxia, embryos 11–15). The heatmap shows relative expression of each transcript (row) and embryo sample (columns). Expression values are count per million (CPM) with clustering by centroid‐linkage. Upregulated genes are in red and downregulated genes are in blue. Examples of genes with differential expression are given in the right‐hand side of the figure. The full list of genes and their level of expression is available as Supporting Information data (Supplement S1)

Transcriptomic analysis

To assess the variation in gene expression across the three treatments, raw counts generated by the Degust analysis were compared between and within treatments. The pattern of gene expression was relatively consistent across the five samples within each of the three treatment groups. In contrast, it differed markedly between the three experimental conditions (0 h, 36 h H and 36 h N). A heatmap of the relative expression of the 250 transcripts that varied most between the conditions showed clearly separated clusters of genes across the three treatments (Figure 1). Clustering of genes was also apparent from a multidimensional scaling plot (Figure 2a). A total of 2422 genes were found to differ (FDR ≤ 0.05, log2FC ≥ ±1) between the 0 h and 36 h N groups, with 1202 upregulated and 1220 downregulated in 36 h N. A total of 3171 genes differed between the 0 h and 36 h H groups, with 1659 upregulated and 1512 downregulated in 36 h H. A total of 4022 genes differed between the 36 h H and 36 h N groups, with 2051 upregulated and 1971 downregulated in 36 h N (Figure 2b, Supporting Information Supplement S1). Overlap between these genes is shown by a Venn diagram (Figure 2b).

FIGURE 2

Differences in gene expression between the three treatment groups. (a) Multidimensional scaling (MDS) plot. Each symbol represents the coordinates (Euclidian distances) for a single embryo collected either immediately after oviposition (0 h, circles), after 36 h incubation under normoxic conditions (36 h N, squares), or after 36 h incubation under hypoxic conditions (36 h H, triangles). Forty‐six percent of the total variance was contributed by MDS dimension 1 while 22% was contributed by MDS dimension 2. (b) Venn diagram showing the numbers of genes with altered expression (>2‐fold difference and p [false discovery rate] < .05) either unique to specific pairwise comparisons of the three treatments or common to two or more comparisons

Gene ontology analysis based on painted turtle orthologues

The ORA analysis was carried out to understand potential functional and pathway level alterations as a consequence of groups of differentially expressed genes. At the broader level of GO parent terms, differences between the 36 h N group and the 0 h group were dominated by groups of genes associated with cellular and multicellular developmental processes and nuclear processes (Figure 3). A similar pattern was seen in the analysis at more specific and fine‐grained level of GO terms (Figure 4). In contrast, differences between the 36 h H group and both the 0 h and 36 h N groups were dominated by groups of genes associated with intracellular organelle structure and function (Figures 3 and 4). Some genes involved in mitosis and cell cycle checkpoints (e.g., sem1, anapc16 and cenpt) were upregulated after incubation in normoxia and downregulated after incubation in hypoxia.

FIGURE 3

Broad‐level, statistically enriched (false discovery rate <0.05), GO terms associated with altered genes in each of the three groups of embryos based on orthologues of the painted turtle. The top five enriched terms for each comparison were combined, and 15 terms are shown for each comparison. BP, biological process; CC, cellular component; GO, gene ontology; MF, molecular function

FIGURE 4

More specific and fine‐grained, statistically enriched (false discovery rate <0.05), GO terms associated with altered genes in each of the three groups of embryos based on orthologues of the painted turtle. The top 5 enriched terms for each comparison were combined, so 15 terms are shown for each comparison. BP, biological process; CC, cellular component; GO, gene ontology; MF, molecular function

Gene set enrichment analysis (GSEA)

A GSEA was performed to understand pathway‐level alterations in the three experimental conditions. Unlike the ORA analysis described above, GSEA utilizes gene level perturbation data to characterize functions or pathways enriched in RNA sequence data. It provides systems‐level context to the list of perturbed genes. Genes involved in DNA replication and mitosis were upregulated in embryos from the 36 h N group relative to embryos at 0 h, while genes involved in glucose metabolism were downregulated, along with genes for specific signalling pathways, such NODAL and adenylate cyclase activation (Figure 5). In contrast, relative to embryos from the 0 h and 36 h N groups, genes involved in mitosis and cell cycle were downregulated in embryos from the 36 h H group.

FIGURE 5

Gene set enrichment analysis. Enriched and depleted HALLMARK pathways for the three pairwise comparisons. (a) Embryos incubated under normoxic conditions for 36 h compared with embryos from freshly laid eggs (36 h N vs. 0 h), (b) Embryos incubated under hypoxic conditions for 36 h compared with embryos from freshly laid eggs (36 h H vs. 0 h) and (c) Embryos incubated under hypoxic conditions for 36 h compared with those incubated for 36 h under normoxic conditions (36 h H vs. 36 h N). Each subplot shows the top 10 up‐ and downregulated gene sets based on the normalized enrichment score (NES). Positive NES indicates that the majority of genes belonging to that pathway were upregulated while negative NES indicates downregulation. Circles show the NES with horizontal lines showing the distance from zero

Expression of genes in the hypoxia inducible factor transcription network

HIFs are master regulators of the biological response to hypoxia (Semenza, 2014) and HIF pathways appear to be important in hypoxia‐induced embryonic arrest in zebrafish (Kajimura et al., 2005, 2006). Thus, we assessed changes in expression of genes under transcriptional control by HIFs. Of the 65 protein coding genes in the HIF‐1α transcription network (Rouillard et al., 2016), 57 were found in the green sea turtle data set. Of the 34 protein coding genes in the HIF‐2α transcription network (Rouillard et al., 2016), 27 were found in the green sea turtle data set. Twelve of these genes are present in both networks. Relative to embryos at 0 h, incubation in normoxia for 36 h led predominately to downregulation of genes in both the HIF‐1 (29 genes, Figure 6) and HIF‐2 (14 genes, Figure 7) networks. In contrast, many HIF‐dependent genes (i.e., HIF‐1 and/or HIF‐2) were upregulated by incubation in hypoxia (Figures 6 and 7). In the HIF‐1α network, 42 of the 57 identified genes were upregulated relative to embryos at 0 h and 37 genes were upregulated relative to embryos incubated for 36 h in normoxia. Similarly, in the HIF‐2α network, 17 of the 27 identified genes were upregulated relative to embryos at 0 h and 16 genes were upregulated relative to embryos incubated for 36 h in normoxia.

FIGURE 6

Heatmap of expression patterns of genes under transcriptional control by hypoxia‐inducible factor‐1. Each row represents a specific hypoxia inducible factor‐1 transcription factor‐associated gene and each column indicates the log2 fold difference for each of the three pairwise comparisons. Blue represents downregulation and the upregulation is indicated by red, with the magnitude of difference shown by the intensity of shading (see inset scale)

FIGURE 7

Heatmap of differential expression of genes under transcriptional control by hypoxia‐inducible factor‐2. Each row represents a specific hypoxia inducible factor‐2 transcription factor‐associated gene and each column indicates the log2 fold difference for each of the three pairwise comparisons. Blue represents downregulation and upregulation is indicated by red, with the magnitude of difference shown by the intensity of shading (see inset scale)

Expression of genes that regulate hypoxia inducible factor transcription activity

We assessed changes in the expression of genes whose products regulate the bioavailability of HIF‐dimers, and thus HIF transcriptional activity. In embryos incubated in hypoxia for 36 h, phd2 (encoding prolyl hydroxylase 2) was downregulated relative to both of the other treatments, vhl (encoding Von Hippel‐Lau protein) was downregulated relative to embryos of freshly laid eggs, while hif1a (encoding HIF‐1α) was downregulated relative to embryos incubated in normoxia for 36 h (Figure 8). In contrast, hif3a (encoding HIF‐3α) was markedly downregulated after 36 h incubation in normoxia relative to the other two treatments (Figure 8). Expression of hif1an (encoding HIF‐1α inhibitor) did not vary significantly between the three treatments (Figure 8).

FIGURE 8

Differential expression of genes that regulate hypoxia‐inducible factor (HIF) activity. Individual black symbols represent each of the five embryos in each of the three treatments groups. Symbols with error bars (red in the online version) represent within‐treatment mean ± standard deviation. Embryos were either removed from the egg soon after oviposition (0 h) or after 36 h of incubation at 26–28°C under normoxic (36 h N) or hypoxic (36 h H) conditions. hif1a encodes HIF‐1α (a), hif3a encodes HIF‐3α (b), vhl encodes Von Hippel–Lindau protein (c), phd2 encodes prolyl hydroxylase 2 (d), and hif1an encodes HIF‐1α inhibitor (e). Pfdr tests for heterogeneity in expression across the three treatment groups after correction for the false discovery rate. When Pfdr ≤ .05, specific contrasts were made using Tukey's test: *p ≤ .05, **p ≤ .01, ***p ≤ .001 for comparison with 0 h; ††† p ≤ .001 for comparison with 36 h N

Expression of genes involved in developmental arrest in other taxa

We assessed changes in the expression of genes identified in previous work to be important in developmental arrest in other taxa. Insulin‐like growth factor binding protein‐1 (IGFBP‐1) mediates hypoxia‐induced embryonic arrest in zebrafish (Danio rerio; Kajimura et al., 2005, 2006). Incubation of green turtle embryos in normoxia was accompanied by a 46% reduction in expression of the igfbp1 gene. In contrast, incubation in hypoxia was accompanied by a 630% increase in igfbp1 expression (Figure 9a).

FIGURE 9

Differential expression of igfbp1 and fam120b. Individual black symbols represent each of the five embryos in each of the three treatments groups. Symbols with error bars (red in the online version) represent within‐treatment mean ± standard deviation. Embryos were either removed from the egg soon after oviposition (0 h) or after 36 h of incubation at 26–28°C under normoxic (36 h N) or hypoxic (36 h H) conditions. igfbp1 encodes insulin‐like growth factor binding protein 1 and fam120b is the turtle homologue of senataxin‐associated nuclease 1 (san1). Pfdr tests for heterogeneity in expression across the three treatment groups after correction for the false discovery rate. When P fdr ≤0.05, specific contrasts were made using Tukey's test: ***p ≤ .001 for comparison with 0 h; ††† p ≤ .001 for comparison with 36 h N

Gene pathways important in developmental arrest in the nematode Caenorhabditis elegans include senataxin‐associated nuclease 1 (san1), encoding SAN‐1, a component of the spindle assembly checkpoint, thus inhibiting progression from the M to G phase of the cell cycle (Miller & Roth, 2009; Nystul et al., 2003). Expression of the homologue of this gene in Chelonia mydas, fam120b (family with sequence similarity 120B; ID: 102946820), was not found to differ significantly among the three experimental groups (Figure 9b)

DISCUSSION

Our results indicate that extending preovipositional embryonic arrest of turtle eggs by incubating them in hypoxic conditions is associated with downregulation of genes involved in DNA replication, cell cycle and mitosis, consistent with continued developmental arrest. It was also associated with upregulation of many HIF target genes, which appear to be predominantly down‐regulated during the first 36 h of normal development after oviposition. Further, the downregulation of genes involved in the cell‐cycle during embryonic arrest in turtles is consistent with forms of developmental arrest in diverse taxa, including yeast (Chan & Roth, 2008), nematodes (Nystul et al., 2003), and teleost fish (Nystul et al., 2003). Thus, even though hypoxia‐induced preovipositional embryonic arrest appears to be unique to turtles among oviparous reptiles, mechanisms underlying hypoxia‐induced developmental arrest per se may be well‐conserved across the diverse taxa in which such mechanisms operate. It is also of interest to note that transcriptional activity in embryos incubated in hypoxia for 36 h differed markedly from that of embryos of freshly laid eggs. This may at least partly explain why hatching success of turtles can be reduced by prolonged extension of arrest by hypoxia after oviposition (Rafferty et al., 2013).

In the current study, the marked differences in transcriptional activity between 0 h and 36 h H groups indicate that the apparent suspension of development in turtle embryos maintained in a hypoxic environment after oviposition is not a static state, but rather an alternative developmental pathway to that occurring under normoxic conditions. Similar observations have been made in the flesh fly (Sarcophaga crassipalpis), in which early and late diapause were found to be transcriptionally distinct (Ragland et al., 2010). The transcriptional differences between the 0 h and 36 h H groups may also at least partly explain why hatching success can be reduced by prolonged extension of arrest after oviposition, and leads to embryonic death if extended beyond a period of a few days (depending on species; Rafferty et al., 2013; Rings et al., 2015; Williamson et al., 2019; Williamson, Evans, Robinson, & Reina, 2017). That is, the extension of embryonic arrest by hypoxia may have adaptive benefits, such as avoidance of adverse environmental conditions during oviposition and synchronization of hatching, at the expense of the viability of developmental processes.

As might be expected, differences in gene expression between embryos from freshly laid eggs and eggs incubated under normoxic conditions for 36 h were dominated by genes involved in normal cellular and multicellular development, with notable upregulation of genes involved in mitosis and the cell cycle. Conversely, genes involved in mitosis and cell cycle were notably downregulated in embryos after 36 h of incubation under hypoxic conditions. This pattern of differences is consistent with the absence of development of freshwater and sea turtle embryos during incubation in hypoxia (Rafferty et al., 2013; Rings et al., 2015; Williamson et al., 2019; Williamson, Evans, & Reina, 2017; Williamson, Evans, Robinson, & Reina, 2017).

HIFs regulate oxygen homeostasis through transcriptional activation (Semenza, 2014), and HIF signalling appears to be a conserved mechanism across a wide range of taxa (Nikinmaa & Rees, 2005; Townley et al., 2017). The range of HIF targets also appears to be relatively well conserved across taxa. For example, in teleosts, as in humans, HIFs regulate expression of vascular endothelial growth factor (VEGF), various glycolytic enzymes, and glucose transporters (Nikinmaa & Rees, 2005), all of which would be expected to contribute to the embryonic response to hypoxia. In mammals, HIF activation by the relatively hypoxic state of the developing embryo plays important roles in placentation and development of the heart and vasculature, lung, and bone (Dunwoodie, 2009). HIF homologues are also important during development in invertebrates, mediating tracheal branching in Drosophila melanogaster and neuronal patterning in Caenorhabditis elegans (Iranon & Miller, 2012). With regard to turtles, hypoxia has been shown to upregulate HIF‐1α in lung epithelial and dermal cells cultured from green sea turtle embryos 3–4 weeks after oviposition (Barlian & Riani, 2018). Mild hypoxia and hypercapnia after breaking of arrest was also shown to slow development in both green sea turtles and loggerhead turtles (Booth et al., 2020). Thus, oxygen‐sensing, including via HIF(s), does appear to operate in the green sea turtle embryo and in the embryos of other turtle species. However, to the best of our knowledge the current work is the first to assess effects of hypoxic incubation to extend embryonic arrest on transcriptional activity in any species of turtle.

Genes that encode proteins in the HIF‐1 or HIF‐2 transcription pathways were predominantly downregulated in embryos incubated in normoxia for 36 h relative to embryos from freshly laid eggs. In contrast, many of these genes were upregulated after incubation in hypoxia for 36 h, even in comparison to embryos from freshly laid eggs. The functional significance of these specific genes in embryonic development in sea turtles has not been established, but our results probably indicate an important role for downregulation of HIF pathways in early post‐ovipositional embryonic development in sea turtles. Whether this includes specific HIF‐dependent genes involved in the maintenance and breaking of preovipositional embryonic arrest remains unknown.

HIF activity is mainly regulated through oxygen‐dependent ubiquitination and proteasomal degradation of the α‐subunits of the HIF dimer (Semenza, 2014). However, changes in the expression of genes encoding HIF proteins (e.g., hif1a and hif3a), the various enzymes involved in degradation of the α‐subunits of the HIF complex (e.g., vhl and phd2), and in genes encoding proteins that modulate HIF activity (e.g., hif1an) could also feasibly influence HIF gene signalling (Semenza, 2014). Incubation of turtle embryos under hypoxic conditions was associated with downregulation of genes that would be expected to both augment (hif1a) and diminish (vhl and phd2) HIF transcriptional activity. These changes may reflect counterregulatory responses to hypoxia‐induced changes in HIF transcriptional activity. Their influence on HIF signalling in turtle embryos and to hypoxia‐induced developmental arrest remains to be determined and merits further investigation.

Hypoxia can induce suspended development in embryonic zebrafish (Danio rerio), at least within 25 h after fertilization, with cells arresting during the S and G2 phases (Padilla & Roth, 2001). The impact of hypoxia on embryonic growth and development in zebrafish appears to be mediated by hypoxia‐induced upregulation of expression of IGFBP‐1 via a HIF‐1 dependent mechanism, resulting in inhibition of IGF‐1‐ and IGF‐2‐stimulated cell proliferation (Kajimura et al., 2005, 2006). Consistent with the plausibility of a similar mechanism operating in turtle embryos, we found that igfbp1 was upregulated, relative to its expression in embryos of freshly laid eggs, when arrest was extended by hypoxic incubation, but not significantly altered after 36 h of development in normoxia.

The transcriptional mechanisms mediating embryonic diapause have also been investigated in the nematode Caenorhabditis elegans (Miller & Roth, 2009). In this model species, embryonic diapause is induced by hypoxia and in turn protects the embryo from its potentially lethal effects (Miller & Roth, 2009). This process appears to be mediated by san‐1, a component of the spindle assembly checkpoint, which inhibits progression from the M to G phase of the cell cycle (Miller & Roth, 2009; Nystul et al., 2003). However, the homologue of the san1 gene in turtles, fam120b, was similarly expressed across the three treatment groups in the current experiment, rendering it unlikely to be a critical factor in hypoxia‐induced developmental arrest in green sea turtles. Nevertheless, in green sea turtle embryos, genes involved in mitosis and cell cycle checkpoints, including sem1, anapc16 and cenpt, were upregulated during normal development in normoxia and downregulated when developmental arrest was extended by incubation in hypoxia. Similarly, in budding yeast (Saccharomyces cerevisiae), genes with roles in ribosome biogenesis, transcription, translation, and many associated processes were significantly downregulated during arrest induced by hypoxia (Chan & Roth, 2008). Thus, while the cellular mechanisms mediating developmental arrest appear to share similarities across diverse taxa, the precise mechanisms mediating this process in the highly conserved and evolutionarily distinct order of testudines (Colston et al., 2020) remain to be determined.

A strength of our study is that our experimental design included biological replicates (n = 5 eggs per group), which appears to be relatively uncommon in molecular ecological research (Todd et al., 2016). Although our sample size was relatively small, we found very clear‐cut differences between the three treatment groups. Importantly, each treatment group included independently generated libraries derived from independent clutches. Thus, we can be confident that our analysis enabled control for variation between clutches and thus maternal identity. However, because female sea turtles are able to store sperm (Pearse & Avise, 2001) and multiple paternity occurs, albeit at a low frequency in green sea turtles at the study location (Fitzsimmons, 1998), we are unable to completely control for paternal identity. Nevertheless, we observed remarkable and consistent differences in gene expression between the three treatment groups, providing confidence in our conclusions.

We acknowledge important limitations of our current study. First, we were unable to build a complete picture of the functional significance of the specific genes and cellular pathways that differed between the three treatments, due largely to the relatively rudimentary functional annotation of the green sea turtle genome. We were also limited by the paucity of available information regarding the roles of oxygen‐sensing mechanisms in embryonic development of turtles and other reptiles. Second, because of the need to freeze embryos as soon as possible after their removal from the egg, we were limited to bulk RNA sequencing rather than newer methods such as single‐cell RNA sequencing. Third, we must also acknowledge that changes in the transcriptome do not necessarily reflect changes in protein synthesis and/or activity (Guppy & Withers, 1999). Thus, our current work must be regarded as hypothesis generating rather than hypothesis testing. Given our current findings, it would also be of great interest, in future studies, to compare the transcriptome of embryos that have recommenced development after extended arrest with that of embryos allowed to develop normally after oviposition. Such a study could provide important insights into the implications, for the developing embryo, of transcriptional processes during extended arrest.

In conclusion, our results indicate that hypoxia‐mediated embryonic arrest in the green sea turtle is accompanied by marked downregulation of genes involved in DNA replication, the cell cycle, and mitosis, as also appears to be the case during developmental arrest in diverse taxa such as yeast (Chan & Roth, 2008) and roundworms (Miller & Roth, 2009; Nystul et al., 2003). It is also associated with marked upregulation of many HIF‐dependent genes, including that for igfbp‐1, which appears to be an important mediator of hypoxia‐induced embryonic arrest in zebrafish. The degree to which igfbp‐1 expression controls the process of preovipositional embryonic arrest in turtles merits further investigation. Our conclusions are limited by the paucity of knowledge regarding the functions of specific genes in embryogenesis in turtles. Nevertheless, to the best our knowledge, the current study provides the first available analysis of the molecular signalling pathways mediating preovipositional embryonic arrest in turtles. This intriguing adaptive strategy, apparently unique to turtle species, has important implications for conservation of endangered turtle species. For example, poor hatching success in the endangered leatherback turtle appears to be mainly due to early embryonic death (Rafferty et al., 2011) and thus may be due to failure of breaking of preovipositional embryonic arrest. Thus, elucidation of the molecular mechanisms mediating embryonic arrest in turtles has the potential to aid conservation strategies.

AUTHOR CONTRIBUTIONS

Richard D. Reina, Angela Gárriz, Roger G. Evans, and Sean A. Williamson conceived and designed the research; Angela Gárriz, Sean A. Williamson, Richard D. Reina, and Roger G. Evans performed the experiments; Angela Gárriz, Anup D. Shah, Deanna S. Deveson Lucas, David R. Powell, and Roger G. Evans analysed the data; Angela Gárriz, Anup D. Shah, Roger G. Evans, Richard D. Reina, Francine Z. Marques, and Sarah L. Walton interpreted the results of experiments; Angela Gárriz, Anup D. Shah, and Roger G. Evans prepared the figures; Angela Gárriz and Roger G. Evans drafted the manuscript. All authors edited and revised the manuscript and approved the final version.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interests.

OPEN RESEARCH BADGES

This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE197628 and https://github.com/MonashBioinformaticsPlatform/Green_sea_turtle_RNA‐Seq.

BENEFIT‐SHARING STATEMENT

Benefits from this research accrue from the sharing of our data and results on public databases as described above.

Table S1

Click here for additional data file.