Sample collection

Total RNA was extracted using RNeasy Kits (Qiagen) [31] from Chrysemys (TSD) and Apalone (GSD). These species were chosen because their sex-determining mechanisms are well characterized [45, 46, 54, 60], and their abundance permits extensive sampling [50]. Embryos were collected at stages 9, 12, 15, 19 and 22 from eggs obtained from a turtle farm and incubated at 26°C and 31°C which are male- (MPT) and female-producing (FPT), respectively, in Chrysemys. Identical incubation conditions were used for Apalone, and these temperatures fall within the optimal thermal range for this species [55]. All procedures were approved by the IACUC committee of Iowa State University. Egg incubation followed previous description [42, 61], using boxes containing moistened sand. Apalone embryos were collected at the same stages as Chrysemys embryos. Embryos reached stages 9, 12, 15, 19, and 22 after approximately 8, 12, 19, 24, and 26 days of incubation at 26°C, and after approximately 5, 9, 15, 19, and 23 days at 31°C, respectively. Total RNA was extracted from trunks (stage 9), adrenal-kidney-gonadal (AKG) complex (stages 12, 15) and gonads alone (stages 19, 22). Stage 15 was targeted as it demarcates the onset of the thermosensitive period before sex determination and differentiation takes place when important differential transcription occurs in the bipotential gonad of mice [62] and of turtle embryos [42]. Stages 9–12 were targeted as important events occur during this time that shape the development of the genital ridge and which are intimately linked to the developing kidneys [63], and because previous studies by us and others show that temperatures experienced prior to the canonical thermosensitive period influence sex ratios and the transcription of genes underlying sexual differentiation in TSD turtles [25, 28, 31, 34, 35, 42, 64, 65]. RNA was quantified using a NanoDrop® ND-1000 Spectrophotometer, and RNA quality was assessed by the presence of ribosomal bands in agarose gels. RNA-seq libraries were generated using pooled samples from ten embryos per temperature per stage for Chrysemys and five for Apalone (a single pool per temperature-by-stage). The resulting 20 libraries (1 library/stage/temperature/species x 5 stages x 2 temperatures x 2 species) were sequenced using Illumina’s HiSeq 2000, and ~35 million 100-bp paired-end reads were obtained per library.

Transcriptome assembly

Reads were splice-mapped across exon boundaries to the Chrysemys genome version 3.0.1 [51] using GSNAP (version 2012-03-23), with the novel-splicing feature turned on [66]. Independently, reads were quality filtered and adapter sequences were removed with Trimmomatic [67] and assembled into species-specific de novo transcriptomes using the Trinity package (release 2013-02-25) [68] and their quality compared with the genome-guided assemblies. De novo transcriptomes were annotated using the Trinotate pipeline [69], mapping the longest open reading frame from each transcript/isoform to the SwissProt protein database [70]. De novo transcripts were mapped to 22,380 genes from the annotated Chrysemys genome [51] using GMAP (version 2012-03-23) [71]. Unannotated transcripts were mapped to Trachemys scripta embryonic transcriptome [72] for comparison using GMAP with default settings. To quantify gene expression levels, the reads from each of the libraries were mapped back to these genes using GSNAP[66]. Then, read-counts for each gene were calculated using HTSeq with the–s (strand-specificity) parameter set to no (which denotes that RNAseq libraries were not strand-specific) [73].

Gene expression normalization

We implemented a novel normalization procedure for read-counts using R version 2.15.2 [74] independently for each turtle species. We employed a mixed approach that combined normalization by the upper-quartile expression levels [75] with normalization to the housekeeping genes Transferrin receptor (Tfr) and hypoxanthine phosphoribosyl transferase 1 (Hprt1) (which were constitutively expressed across all stages in both species). This approach permitted validation of the transcriptomic expression levels by comparison to extensive expression data from Chrysemys obtained by qRT-PCR of candidate genes from individual embryos, and which were normalized to housekeeping genes during a separate study [42] using independent RNA samples from the present study. To test the effect of the normalization procedure on the number of differentially-expressed genes, we conducted multiple Fisher exact tests between transcript expression levels that were normalized by (1) upper-quartile only (procedure 1, named UQ100), (2) upper-quartile after eliminating the top 1 percentile of transcripts with the highest expression (procedure 2—UQ99), (3) upper-quartile and house-keeping gene normalization (procedure 3—UQHK100), and (4) upper-quartile and house-keeping gene normalization after eliminating the top 1 percentile of transcripts with highest expression (procedure 4—UKHK99).

Differential expression tests

Differential expression tests were performed per developmental stage between the MPT and FPT for Chrysemys (TSD), which correspond to high/low temperature for Apalone (GSD). Fisher’s exact tests were used to evaluate differential gene expression between high and low temperature treatments at each developmental stage in each species, as this statistical method permits the analysis of un-replicated samples [75–77]. Briefly, Fisher’s exact test is based on a 2x2 contingency table (Table 2)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Contingency table for Fisher’s exact test.

https://doi.org/10.1371/journal.pone.0172044.t002

where n_ki is the observed read count for focal gene X (k = 1) or for all other genes in the transcriptome (k = 2) for treatment i (i = 1 for 26°C and i = 2 for 31°C); n₁₁ + n₁₂ is the marginal total for focal gene X; n₂₁ + n₂₂ is the marginal total for the remaining genes in the transcriptome; n₁₁ + n₂₁ is the marginal total for the 26°C treatment; n₂₁ + n₂₂ is the marginal total for the 31°C treatment; N is the grand total (reviewed in [76]). This approach tests the null hypothesis that H₀: θ = 1, where θ = π₁₁π_22/π₁₂π₂₁ and where π_ki is the true proportion of counts in cell k,i (k = 1 for focal gene X or 2 for all other genes; i = 1 for 26°C and 2 for 31°C). In other words, this approach tests the null hypothesis that the proportion of counts (gene expression) of focal gene X between 26°C and 31°C is the same as for all other genes. It is noted that the lack of replication implies that differences between treatments detected here are not generalizable in the way that inferences from replicated datasets would be, and thus, they should be treated as suggestive hypotheses that warrant corroboration by further analyses. The resulting p-values of Fisher’s exact tests were corrected for false discovery (FDR) [78]. Then, we concentrated on the highly differentially expressed genes meeting a stringent FDR-corrected cutoff of 1e^-10 chosen arbitrarily. Differentially-expressed genes were annotated using the KEGG database [79]. Enrichment analyses of Gene Ontology (GO) categories against their respective species-specific transcriptomes were conducted using the DAVID Bioinformatics knowledgebase [80]. Additionally, turtle transcriptomes were tested for enrichment against the mouse gonadal transcriptomes [21] using DAVID [79] to test for transcriptional divergence between the turtle and mammalian lineages. In a complementary approach, we randomly subdivided each read library into 2 and 3 subsets (or “subsamples”) [81], identified the GSNAP alignments corresponding to these subsamples, and regenerated read counts per gene. We then used DESeq [82] and EdgeR [83] to independently determine the differentially-expressed genes by leveraging the multiple subsamples while controlling false discoveries at 1%. Finally, the R package WGCNA [84] was used to identify modules of genes co-expressed across turtle embryonic stages in the original set of libraries as well as the subsamples.

Transcriptome assembly

RNAseq data was obtained from Chrysemys embryonic tissues before the thermosensitive period when the gonads are not yet discernible (stages 9 trunks, stage 12 adrenal-kidney-gonad complexes—AKGs), at the onset of the TSP when gonads are hard to separate from surrounding tissue (stage 15 AKGs), and during the mid and late TSP (stage 19 and 22 gonads), from male-producing temperature (MPT = 26°C), and female-producing temperature (FPT = 31°C) [42]. Identical incubation conditions and sampling scheme were followed for Apalone. Such mix-tissue sampling has been used in recent transcriptomic studies of other TSD turtles [25]. De novo transcriptome assemblies constructed using Trinity [68] resulted in a high percentage of mapped reads (>92%), and high representation of Core Eukaryotic Genes (CEGs) (>77%) and mammalian urogenital development pathway genes in both species (>96%) (Table 3). All subsequent analyses reported here are based on the de novo transcriptome assemblies. We also tested the alternative approach using reference genome-guided assemblies. However this approach was discarded because while the mapping rate of the Chrysemys libraries to the Chrysemys reference genome [51] was high (97%), that of Apalone was only 44% (Table 3) and resulted in significantly fewer gene models for Apalone (5,596 unique to Apalone, 14,661 unique to Chrysemys, and 23,465 overlapping). The problem persisted when using as a reference the genome of Pelodiscus sinensis [85], a close relative of Apalone, because the P. sinensis assembled genome is more fragmentary than that of Chrysemys, as evidenced by the lack of complete exonic sequences for several genes (such as some homologs of mammalian urogenital genes: Ctnbb1—missing part of exon 4; Wnt4—missing exon 1 and part of exon 4; Dmrt1—exon 1 is mis-assembled/mis-annotated; Sox9—exon 1 is mis-assembled).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Genome-guided and de novo transcriptome assembly results for Chrysemys picta (CPI) and Apalone spinifera (ASP).

CEGs = Core Eukaryotic Genes.

https://doi.org/10.1371/journal.pone.0172044.t003

Normalization schemes to identify differentially-expressed genes per species

Gene read abundance was normalized multiple ways, first to the housekeeping genes Tfr and Hprt1 which were constitutively expressed across all developmental stages in both Chrysemys and Apalone, and then by the standard upper-quartile normalization [75]. The order of these steps did not affect the overall assessment of gene expression. When compared to other normalization procedures described in the methods (UQ100 = upper-quartile only, UQ99 = upper-quartile after eliminating the top 1 percentile of transcripts with the highest expression, UQHK99 = upper-quartile and house-keeping gene normalization after eliminating the top 1 percentile of transcripts with highest expression), the chosen UQHK100 approach resulted in fewer differentially-expressed genes than using the upper-quartile alone (Fig 1A), and therefore, it is a more conservative approach. Furthermore, UQHK100 normalization revealed differential expression patterns which were most consistent with extensive qPCR data of several candidate genes previously obtained for Chrysemys using a completely independent set of RNA samples [42], as determined qualitatively by visual inspection of the expression profiles over developmental time for individual genes to identify broadly concordant expression patterns (examining whether differential expression was present and in the same direction per stage and gene between the two studies) (Fig 2). Therefore, we used UQHK100 to identify differentially-expressed genes for further enrichment analyses to ensure unbiased comparisons between species.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

Overlap of differentially-expressed genes by temperature (controlling FDR at 1e-10) for various categories: (a) overlap between male- and female-producing temperatures at Stage 22 of Chrysemys based on expression levels normalized by (1) upper-quartile (UQ100); (2) upper-quartile excluding the top 1 percentile of transcripts with the highest expression (UQ99); (3) upper-quartile and house-keeping genes (UQHK100); and (4) upper-quartile and house-keeping genes excluding the top 1 percentile of transcripts with highest expression (UQHK99). (b) and (c): across stages in Chrysemys picta and Apalone spinifera. (d): Overlap of genes of interest present in turtles with the mouse gonadal genes described in Jameson et al., (2012). DE = differentially-expressed, GoI = genes of interest described in Table 5. Differential expression was assessed between temperatures at each developmental stages for each species separately.

https://doi.org/10.1371/journal.pone.0172044.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Average qPCR expression at 26°C (blue) and 31°C (red) across developmental stages from qPCR experiments (top panels; modified from Valenzuela et al., 2013) and RNAseq (bottom panels; this study).

Mean and one standard deviation are presented. Asterisks denote significant differential expression by temperature. Boxed stages fall within the thermosensitive period.

https://doi.org/10.1371/journal.pone.0172044.g002

Gene annotation and genes of interest

The longest open-reading frame was used to predict protein sequences, and between 26% and 28% of the transcripts from the de novo transcriptome assembly per species were represented in the SwissProt protein database. The Chrysemys transcriptome showed an overall higher representation of annotated genes in the Chrysemys genome [86] than the Apalone transcriptome (Table 4). Of the transcripts not annotated in SwissProt, 252 from Chrysemys and 169 from Apalone correspond to non-coding RNA sequences (ncRNAs) [87] as identified by BLAST [88]. Among these, 115 out of the 252 transcripts In Chrysemys and 70 out of the 169 in Apalone were differentially-expressed (S2 Table). A majority of these ncRNA transcripts were upregulated at 26°C in both species of turtles (Chrysemys’ MPT), 36 were upregulated during developmental stages 9, 19 and 22 in Chrysemys and another 38 were upregulated only during stages 19 and 22 (Chrysemys’ TSP). Of the 25 ncRNA transcripts that overlap, all but 3 show similar differential expression pattern in both species (S2 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Number of annotated genes identified in the transcriptomes of Chrysemys picta (CPI) and Apalone spinifera (ASP) by treatment (26°C and 31°C) and developmental stage (9, 12, 15, 19, 22) out of the 22,380 annotated genes present the Chrysemys picta reference genome.

https://doi.org/10.1371/journal.pone.0172044.t004

To search for candidates which may be involved in the conversion of temperature signals into sex-specific development and thus may have a significant role in TSD based on their differential expression pattern by temperature, we focused on known genes involved in (a) vertebrate sex determination/differentiation in model mammals and birds [18, 42, 56, 89], (b) epigenetic modification [57], (c) hormonal pathways [58] and (d) general sensing responses [59], out of the variety of annotated genes. These target categories included heat-shock genes, transient receptor potential genes, germ-line and histone-related genes, androgen- and estrogen related genes and genes linked to human/chicken sex chromosomes (http://www.ensembl.org/info/data/ftp/index.html). While the overall composition of the transcriptomes of the two turtle species was similar with regard to these categories, a few genes present in Chrysemys’ transcriptome were notably absent in Apalone across all stages, including some genes X-linked in human and Z-linked in chicken (Table 5). Further, some genes were differentially expressed by temperature in both turtles (S2–S11 Tables). Interestingly, a number of genes that are involved in histone modification show low temperature bias (MPT) just before the onset of the thermosensitive period in Chrysemys but high-temperature bias in Apalone: histone H1-x like protein (H1x-like), histone chaperone protein (Asf1B-like), H3-Histone family 3B (H3f3b) and Nuclear autoantigenic sperm protein (Nasp). Details about the transcriptional response of these genes of interest are presented in the discussion.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Number of genes per category of interest present in the transcriptomes of Chrysemys picta (CPI) and Apalone spinifera (ASP).

https://doi.org/10.1371/journal.pone.0172044.t005

Differential expression in painted turtles (TSD)

The differential expression patterns of five candidate genes previously detected by qPCR [42] and used here for validation of our transcriptomes, were recapitulated by the RNA-seq data when qPCR differences in expression were large. However, smaller expression differences identified by qPCR at certain embryonic stages passed undetected in our RNA-seq transcriptomes (Fig 2). We then chose a highly conservative p-value cutoff value of 1e^-10 to further correct for false positives in our differential expression tests. This approach revealed significant overlap of highly differentially-expressed genes across developmental stages for both species (Fig 1B and 1C). The same differentially-expressed genes were also recovered with a second approach where reads from each RNAseq library were randomly subdivided into multiple representative subsamples [81] which were then used in the differential expression analysis using both DESeq and EdgeR toolkit [82, 90]. Among these, we identified 1065 genes that were differentially-expressed only in Chrysemys across development (S3 Table). Some results of particular interest are highlighted below.

As mentioned above, our Chrysemys RNA-seq data corroborated qPCR results for five known gene homologs involved in gonadogenesis in Chrysemys, namely Sf1, Wt1, Sox9, Aromatase and Dax1 [42] which serve to validate our transcriptomes (Fig 2). This approach is similar to the validation level of a recent turtle transcriptome study [91]. RNA-seq was also used to profile a large number of candidate sex-determining genes in Chrysemys and Apalone (Figs 3 and S2), which were previously uncharacterized in turtles, contributing extensively to our understanding of the full composition of the turtle sex-determining regulatory network. A significant number of these genes in Chrysemys show MPT-bias late in the TSP (Stage 22), two (Igf1r and Insr) before the formation of the bipotential gonad (Stage 9), and only one (β -catenin—Ctnnb1) shows FPT-bias before the thermosensitive period (Fig 4). Overall, the known vertebrate sex-determination/differentiation network exhibited greater responsiveness to temperature in Chrysemys with expression profiles during the thermosensitive period being more concordant with those predicted by the function of these genes in mammals (i.e., genes linked to testicular formation show higher transcription at MPT and those linked to ovarian formation in mammals show higher transcription at FPT in Chrysemys, except for Esr1 and Hsp90ab1 [Figs 3 and 5]). Contrastingly, expression patterns in Apalone were less thermosensitive and more variable, shifting between MPT-bias and FPT-bias throughout development [Figs 3 and 5]. We note that gonads alone were not available from the earlier-stage embryos (we used stage 9 trunks and stage 12/15 AKGs), either because they have not discernibly developed yet (stages 9, 12), or because they cannot be confidently separated without contamination from the surrounding tissue (stage 15). Nonetheless, tissues from these early embryos (which do contain the developing gonad) represent a critical sampling time to detect early responses of key candidate genes to temperature, and their transcriptional profiles revealed that the machinery underlying sexual development is active well before the “canonical” temperature-sensitive period.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Differential expression in Chrysemys picta and Apalone spinifera turtles based on RNAseq of a subset of genes involved in the mammalian urogenital pathway and other genes of interest for turtle gonadogenesis.

The list of genes and activation or repression information derives from Valenzuela 2008b, Liu et al, 2009, Chassot et al., 2012, Eggers et al., 2014, Lai et al., 2014 plus other sources cited in the text, while information about the timing and pattern of expression correspond to those observed in turtles during the present study. Approximate equivalency is provided between mice and turtle developmental stages of gonadal development. [not sig. diff. exp. = Not significant differential expression].

https://doi.org/10.1371/journal.pone.0172044.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. P-values of differentially-expressed genes (after applying Benjamini-Hochberg correction) linked to mammalian urogenital pathways, showing higher expression at 26°C (blue) and 31°C (red).

Highly differentially-expressed genes (while controlling false discoveries at 1e-10) identified in dark blue and dark red. Light colored cells denote significance at a standard α = 0.05. GA = gene absent in transcriptome.

https://doi.org/10.1371/journal.pone.0172044.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Eigengene networks and dendrograms illustrating the co-expression patterns in developing gonad of (a) Chrysemys at 26°C (b) Chrysemys at 31°C, (c) Apalone at 26°C and (d) Apalone at 31°C, of 981 genes of interest described in Table 5. Colors along the x and y-axes represent clusters of genes showing similar expression. Gene order varies by plot along the X and Y-axes.

https://doi.org/10.1371/journal.pone.0172044.g005

Gene Ontology (GO) Enrichment

At each embryonic stage the differentially-expressed genes were enriched for a number of GO pathways. While no GO pathways were consistently deployed in both turtle species throughout development, GO pathways relating to translation and translational elongation were present across three out of the five stages in both taxa. Overall, we found more shared pathways by stage (including general cell functions such as mitotic cycles, mRNA processing and RNA-splicing) between Chrysemys and Apalone before stage 15 (onset of TSP in Chrysemys) than later in development, suggesting that temperature triggers different network modules after the onset of thermosensitive period in Chrysemys than in Apalone. Indeed, enriched pathways during stages 9 and 12 in both turtles, including intracellular transport, protein localization and protein catabolic processes remained enriched only in Chrysemys after stage 15 (S4 Table). In contrast, genes upregulating protein ubiquitination and ubiquitin protein ligase were enriched only in stage 12 of Apalone and not in Chrysemys.

Novel transcripts

Around half of Chrysemys (53%, or 150,195/279,903) and Apalone transcripts (54% or 152,579/279,753) were absent in SwissProt or ncRNA databases. However, 87% (131,131) of Chrysemys transcripts were mapped to the Chrysemys genome, of which only 7% (10,660) were unannotated, indicating a gap between the SwissProt/ncRNA databases and the annotated Chrysemys genome. Among these, most differentially-expressed transcripts are MPT-upregulated at stages 9, 19 and 22 (S1 Fig). Only 50% of Apalone novel transcripts could be mapped to the Chrysemys genome and 26% to the Trachemys scripta whole-body transcriptome [72], while 68% of Chrysemys novel transcripts mapped to T. scripta. Since this mapping was carried out under low stringency conditions, this difference is likely due to the absence of many novel transcripts in Apalone, the GSD turtle, as extensive divergence is not expected in coding sequences between Chrysemys and Apalone given that transcripts from Chrysemys and Apalone are more often identical than not (data not shown).

Gene clustering and co-expression

The genes of interest (Table 5) clustered into modules of co-expression patterns across embryonic stages that differ between turtles. Interestingly, Chrysemys showed stronger clustering differences between temperatures, with more gene modules at 26°C (45) than at 31°C (10) (Figs 5A, 5B, S3A and S3B). In contrast, Apalone did not differ significantly between temperatures in the number of gene modules (Figs 5C, 5D, S3C and S3D). Similar species-specific clustering differences were also detected for 157 core eukaryotic genes (data not shown). The composition of the largest clusters also varied by temperature in terms of GO pathways (Table 6).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Categories of GO pathways enriched (at p = 0.05) in the largest co-expressed clusters in the Chrysemys and Apalone embryonic transcriptomes.

https://doi.org/10.1371/journal.pone.0172044.t006

Transcriptome assembly

The genome-guided transcriptome assembly using Chrysemys as reference [51] worked well for Chrysemys but produced poor results for Apalone (44% mapped reads and fewer gene models; Table 3), underscoring the extensive divergence accrued in these turtle genomes since their lineages split >180mya [12]. Using the Pelodiscus sinensis genome as reference [85] did not solve this problem, and comparative approaches require common analyses for all species. However, de novo transcriptome assemblies had similar high quality and permitted the discovery of novel transcripts previously unidentified in public databases.

Contrastingly, the Chrysemys genome [51] was useful to map the transcripts from the de novo assemblies (since transcripts are longer and align better than reads) to quantify the representation of annotated genes per library (Table 4). The Chrysemys transcriptome had an unsurprising slightly higher representation of annotated genes overall (62%) than Apalone (57%). The P. sinensis genome was excluded here because its annotation is less extensive than Chrysemys’. These observations highlight the need to improve current turtle genome assemblies [52], whose re-annotation is aided by RNAseq data, and to sequence additional genomes from representative phylogenetic lineages to illuminate turtle and vertebrate genome evolution and aid future ecological and evolutionary genomic studies.

Using these assembled transcriptomes we address several hypotheses about (1) the composition of the vertebrate gene network governing sexual development in turtles, (2) how it might perceive and transmit environmental temperature inputs, (3) its similarities or differences between TSD and GSD, and (4) whether it may contain reptile-specific elements, as described in the sections below.

1. The vertebrate gene regulatory network underlying sex determination/differentiation is present and active in TSD and GSD turtles.

The following paragraphs highlight the transcriptional patterns of known vertebrate determination/differentiation genes found in our transcriptomes (full gene names are presented in Table 1).

Genes in the vertebrate gonadal network but unknown in turtles

RNAseq provided novel embryonic transcriptional profiles in Chrysemys and Apalone of several vertebrate genes unstudied in painted and GSD turtles (Fig 4), including genes implicated in testicular differentiation in mammals [Amh, AR, Cbln4, Dhh, Dmrt2, Fgf9, Fgfr2, Fhl2, Fog2, Pgds, Ptch1, Srd5a2 and Vnn1], genes involved in ovarian formation [Ctnnb1, Esr2, Foxl2, Gata2, Rspo1 and Wnt4], and genes important for both testicular and ovarian function or general gonadogenesis prior to their sexual commitment [Cbx2, Ck1, Gsk3b, Apc, Insr, Igf1r, Kdm3a, Six1, Six4, Dmrt3, Emx2, Esr1, Gata4, Lhx1, Lhx9] [31, 89].

Known genes in the turtle gonadal network

Multiple genes involved in turtle sex determination have been studied with candidate gene approaches, including Vasa, Dazl, Amh, Foxl2, Dmrt1, Aromatase, Androgen receptor and Estrogen receptors α and β, among others [38, 94–97]. Our transcriptomes recapitulated expression patterns from qPCR for five such genes previously reported in Chrysemys [27, 28, 31, 34, 35, 98], although subtle differences passed undetected in the transcriptomes (Fig 3), perhaps because transcriptomic inferences have lower power overall than qPCR approaches [43, 44]. We note that to avoid false positive results from the absence of biologically replicated transcriptomes, and from the potential bias introduced by the lower number of embryos per RNA library for Apalone compared to Chrysemys in our study, we applied a stringent cutoff of 1e-10 to control for false discoveries, and discarded any genes with lower significant differential expression (Fig 5). Taken together, our results highlighted below reveal that the vertebrate gene network regulating primary sexual development is highly conserved in its composition and is active in turtles, but regulated differently between TSD and GSD turtles, and between turtles, mammals and birds (Fig 2, S1 Table).

Importantly, we found that overall, Chrysemys exhibits more extensive thermosensitive transcription of the genes in the vertebrate regulatory network than Apalone. Also, the expression profiles during the thermosensitive period are fairly concordant with those predicted by the function of these genes in mammals, such that numerous elements involved in mammalian testiculogenesis show higher transcription at MPT and elements involved in mammalian ovariogenesis show higher transcription at FPT in Chrysemys (Fig 5). The notable exceptions are Esr1 and Hsp90ab1, which display a transcription pattern opposite of that in mammals, perhaps revealing a new function of these genes in turtles. Such turnovers are not unprecedented as extensive transcriptional evolution has taken place in elements of this network across vertebrates [42]. In contrast, expression patterns in Apalone were more variable, shifting between MPT-bias and FPT-bias throughout development, as might be expected under relaxed selection after the evolution of GSD from TSD.

2. How are thermal signals converted to TSD sex-specific development? Does thermal-signal transmission involve the epigenetic machinery, hormonal pathways, and general sensing responses?

Differentially-expressed genes by temperature

We searched for genes showing thermosensitive expression in Chrysemys (TSD) to uncover candidate temperature sensors or transducers that may activate TSD male and female gonadogenesis. We detected many such genes, including numerous homologs of mammalian sex determination/differentiation genes, plus previously undescribed candidates (Fig 5) as mentioned above. A comparison of differential expression by temperature across select vertebrates is summarized in S1 Table. It is important to note that trunks were examined at stage 9, AKGs at stages 12 and 15, and gonads alone at stages 19, 22, such that gene expression from the developing AK or other organs may contribute to the differences or lack thereof between temperatures at these earlier stages of development. Of the genes in this regulatory network, Amh, Ar, Esr1, Fog2, Gata4 and Lhx9 show significant MPT (26°C) bias at stage 22 in Chrysemys (i.e., during the thermosensitive period when temperature has a strong effect to bias the resulting sex ratios) and are thus important candidate genes for turtle thermosensitive testicular differentiation in TSD vertebrates that deserve further functional research. Noteworthy, our analyses show a consistent FPT (31°C) bias in the expression of the Cold-inducible RNA-binding protein (Cirbp) in both Chrysemys and Apalone. This pattern of temperature-dependent expression in the Cirbp gene was shown to help govern ovarian development in the snapping turtle Chelydra serpentina (a TSD turtle) via allelic dimorphism, and represents an important novel TSD candidate gene [40]. Combined, our results and those in Chelydra indicate that Cirbp upregulation at warm temperature might be ancestral to turtles and thus perhaps relict in Apalone (as proposed for Wt1 and Dax1 [27, 31]). Alternatively, the differential expression of Cirbp in Apalone, identical to that of Chrysemys and Chelydra, would call into question its purported key role as a TSD element [40]. We note that in the absence of sex information from the Apalone embryos in our study, it remains plausible that some of the thermal effects detected here might be due to sampling effects (e.g., to the predominance of one or the other sex in some of the pools of RNA from Apalone used in this study). Further research is warranted using recently developed sexing techniques for Apalone [93] to test these alternative hypotheses directly.

Genes in other functional categories

We explored additional functional gene categories of plausible transducers of the temperature signal to gonadal developmental, some previously known as temperature-sensitive or linked to gonadal formation in other animals [57–59]. These include vertebrate genes involved in gonadal and germ-line differentiation, androgen- and estrogen related genes, and genes linked to sex chromosomes, heat-shock and transient receptor potential genes and histone-related genes. Many of these genes exhibited thermosensitive expression in both turtles (S2–S11 Tables), including genes involved in histone modification, several kinases, genes involved in androgen- and estrogen signaling pathways, sex-linked genes and heat shock proteins. Overall, transcriptome composition was similar between species with some noticeable differences. Namely, Apalone’s transcriptome exhibited slightly lower representation of kinases, ubiquitin- and histone-related genes (Table 5), although we confirmed that they are present in the genome of Apalone using BLAST. Kinases are indispensable for cell functioning and orchestrate many cellular processes. One of these, the protein kinase Map3k4, directly affects Sry and Sox9 expression in bipotential mice gonads, inducing testicular development [99, 100]. Several heat shock proteins show sexually-dimorphic expression in American alligator (TSD), potentially influencing sex determination [59], but only one showed the same pattern in alligator and Chrysemys (Hspb1—FPT-biased) whereas the heat shock gene Hsph1 was upregulated at MPT in the alligator and at FPT in Chrysemys (stage 12) (S6 Table). Also surprisingly, comparisons of our turtle and the recent alligator transcriptomes [24] revealed only a handful of genes (S6 and S7 Tables) with shared MPT- and FPT-specific expression pattern between alligator and Chrysemys [only 31 out of 207 genes reported to be upregulated at FPT, and 43 out of 250 reported to be upregulated at MPT in the alligator were shared with the Chrysemys gonadal transcriptome (S6 and S7 Tables)]. The expression pattern of a few genes was similar between the alligator and Chrysemys either by sex or by relative temperature. For instance, Sf1, Amh and Amhr2 were upregulated at MPT during the TSP in both species [stages 19 (Sf1, Amh) and 22 (Sf1, Amh and Amhr2) in Chrysemys], while Crabp2 and Hspb1 were upregulated at FPT in both taxa. Finally, differences in expression were observed in the transient receptor potential gene Trpc4ap and Tex15, which were both upregulated at MPT in the alligator but not in Chrysemys. No other genes in our categories of interest discussed above overlapped between the Chrysemys and alligator. These observations support the notion that significant transcriptional evolution has accrued even among TSD regulatory networks [42] and leads to the hypothesis that perhaps numerous temperature-specific responses may be species-specific.

Sex-linked genes such as Nf2 and Prdx4 [101, 102] are differentially-expressed by temperature in Chrysemys at stage 9. Similarly, we found thermosensitive expression for Serpinh1, Hsp90ab1 and Hspa8 across stages in both turtles (S8 Table) (and potentially relic in Apalone) while expression is monomorphic in mouse [21], suggesting their potential turtle-specific role in gonadogenesis. Additional genes differentially-expressed in turtles but at relatively later stages in the mouse gonad [21] include Ctnnb1, a known ovarian inducer [89] (which was early acting at 31°C in both Chrysemys and Apalone) and Git2 (a sex-linked gene in Pelodiscus sinensis [101]), among others (Fig 1D, S3 and S16 Tables). This indicates that genes already important for gonadal formation at downstream stages in mammals, such as Ctnnb1, have been recruited for an earlier than anticipated temperature-specific response in turtles, underscoring the extensive ontogenetic evolution that has accrued in the transcriptional patterns of this regulatory network in vertebrates in general and in turtles in particular [25, 42].

We also identified numerous differentially-expressed ncRNAs (S2 Table) by temperature, which interestingly, exhibited upregulation at the colder temperature (26°C) in the majority of cases and in both turtles, although few ncRNAs transcripts overlapped between species. The biological function of these ncRNAs remains an open question worthy of future study.

Thermosensitive response of signaling pathways

Distinct cell types may derive from a handful of cell signaling pathways [103]. We found evidence that numerous signaling pathways are differentially regulated by temperature. For instance, Jak-Stat signaling, involved in cell proliferation and hematopoiesis [104] exhibits male-bias of Egfr expression at stage 22 in Chrysemys. Nf-κB signaling plays a role in immune and stress response [105], and involves members of the hypoxia-induced Tumour necrosis family (Tnf) [106] and Breakpoint cluster region (Bcr) [107] which showed female-bias in Chrysemys at stage 22 and male-bias at stage 9, respectively. The receptor gene Vegf, which regulates sex-specific gonadal vasculogenesis [108] was also female-biased at stage 15 in Chrysemys and high-temperature biased during stages 19–22 in Apalone. Further, retinoic acid has been identified to induce meiosis in mice germ cells regulating ovarian formation [109]. Two retinoic acid binding proteins were differentially transcribed: Crabp1 (male-biased during Chrysemys TSP) and Crabp2 (female-biased pre-TSP and during the TSP in Chrysemys, and upregulated at high-temperature from stage 15 onwards in Apalone). Among the signaling pathways implicated in vertebrate sex determination, Foxl2 and members of the Wnt signaling pathway regulate ovarian formation [110, 111]. Wnt activates Ctnnb1, which inhibits Sf1 from activating Sox9 and inducing testiculogenesis [112]. The canonical Wnt machinery including Ck1, Apc, and Gsk3 show male-bias during stage 9 in Chrysemys, and monomorphic expression in Apalone, indicating that Wnt signaling is active in TSD and GSD turtles, but deployed differentially by temperature. Members of the Mapk signaling family, required for Sry activation testiculogenesis in mice [100] were also low-temperature biased in turtle bipotential gonads despite the absence of Sry (Map3k3 at stage 9 in Chrysemys; Map3k7 at stage 12 in Apalone) (S9 Table), rendering them additional candidates for functional tests. Akt signaling is directly activated by Fgf9 in mice, promoting steroidogenesis (Lai et al. 2014), but neither gene showed thermosensitive expression. Finally, Ras-mediated signaling is implicated in sex myoblast migration in nematodes [113], and a subtle thermosensitive expression was detected in Chrysemys (S9 Table). Thus, Jak-Stat, Nf-κB, retinoic acid, Wnt, and Mapk signaling are potentially involved in TSD gonadogenesis, while this process appears independent of Akt and Ras-mediated signaling.

3. What makes gonadal networks sensitive to environmental temperature in TSD and not in GSD?

Gene enrichment analysis

Enrichment analyses of Gene Ontology (GO) categories represented in the transcriptomes using DAVID [80] revealed that the two turtle species examined here shared more pathways before stage 15 overall, except for chromatin organization and chromatin modification pathways which were enriched only at stage 9 in Chrysemys, and pathways linked to protein ubiquitination and ubiquitin protein ligases which were enriched only at stage 12 in Apalone. Ubiquitination is a post-translational modification that results in protein degradation [114]. This suggests that temperature triggers a different set of downstream pathways in Chrysemys than in Apalone potentially leading to sexual fate determination by temperature in the former. Pathways including intracellular transport, protein localization and protein catabolic processes, were enriched in both species at stages 9 and 12, but remain enriched only in Chrysemys after stage 15 (S4 Table). Thus, important evolutionary changes may have occurred in GSD turtles in the machinery underlying gonadogenesis preceding the thermosensitive period of TSD turtles (before stage 15), perhaps inactivating genes regulating the male- and female-specific TSD pathways and consequently, determining sex independent of temperature. These and earlier findings (reviewed in [42]) underscore that key thermosensitive events for sexual development occur in early embryogenesis and deserve further research.

Finding differential expression prior to the onset of the canonical thermosensitive period in Chrysemys is of particular importance as any such gene may be the key TSD master element that senses the environmental temperature signal or a key activator of the thermosensitive period. Additionally, the disruption of these potential TSD master elements could lead to GSD evolution. Notably, Ctnnb1 (β-catenin, a member of the Wnt signaling pathway) showed female-bias at stages 9 and 12 in Chrysemys consistent with its involvement in early ovarian formation in mammals [115, 116]. Similarly, Ctnnb1 shows high-temperature bias in Apalone during stages 12 and 15 but not at stage 9. Follistatin (Fst) a gene activated by Ctnnb1 in mice bipotential gonads (Eggers et al. 2014) showed slight female-bias (α = 0.05) in Chrysemys at stage 9, and slight high-temperature bias in Apalone at stage 15 (α = 0.05), suggesting that Ctnnb1 could also activate Fst in turtles. Further, this results indicates that Ctnnb1 and Fst thermosensitive expression may be ancestral to cryptodiran turtles (the suborder to which Chrysemys and Apalone belong), which if true, would indicate relic thermosensitivity in Apalone, and underscores that downstream elements may be key to rendering GSD immune to temperature as is the case of Wt1 and Dax1 whose thermoresponsive activity is rendered ineffective via the loss of thermosensitivity of Sf1 [27, 31]. Indeed, Ctnnb1 is a repressor of Sf1 [112], a gene with thermosensitive expression in Chrysemys but not Apalone [28]. Also noteworthy, Chrysemys (and not Apalone) shows upregulation at MPT of Insr and Igf1r during stage 9, two genes indispensable for testiculogenesis in mice [117] which antagonize the Ctnnb1-Wnt signaling pathway essential for ovarian formation. This indicates that the same molecular antagonism exists in TSD turtles, is active before the canonical thermosensitive-period, and could influence growth trajectories via the insulin receptor family, inducing male determination [118] and other sexual dimorphisms with potential temperature-specific fitness consequences [119]. These early differences between TSD and GSD systems may have functional significance for the evolution of sex determination.

Comparison with the slider turtle Trachemys scripta

Both Trachemys and Chrysemys have TSD and diverged near 30 million years ago [120]. qPCR evidence exist that transcriptional patterns of some of genes involved in sexual development have diverged between these two species during this time [42], but the extent of differences and similarities in other elements of this network remains unclear. Here we compared the expression stage-by-stage of key common elements between our transcriptomes and the recent Trachemys transcriptome [25] (Fig 6). We observed similarities in the expression of certain genes between species, in terms of upregulation at MPT or FPT in both the embryo and gonad simultaneously, or in the gonad alone (Fig 6). It should be noted that the Trachemys study did not sample stage 9 as in our study, whereas stages 17 and 18 were examined in Trachemys and not in our study. Additionally, the thermosensitive period differs slightly between Trachemys (stages 15–21) [121] and Chrysemys (stages 16–22) [46], such that while numerically different, both studies contain data from stages at the onset and at the end of the thermosensitive period that can be compared (Fig 6). Differential transcription of numerous genes in this subnetwork was observed in Trachemys at stages 12 and 15 that were not observed in Chrysemys either because the timing of their expression is divergent, or because they passed undetected by our unreplicated transcriptomic approach, yet when expression was detected the pattern was concordant between species except for Avil (Fig 6). Greater similarities were observed later in development in both the embryo and gonad during stage 19 at FPT (Anpep, Cutc, Des, Dgka and Eif4a2). Among genes that were not differentially expressed during stage 12 but showed MPT or FPT upregulation at stage 15 or later (and thus are candidates for a role as temperature-triggered genes), similar profiles were observed late in the TSP for Aromatase, Avil, Dpt, Mertk, and Twist1, all of which were upregulated at FPT during stage 21 in Trachemys and stage 22 in Chrysemys. Genes upregulated at MPT in both turtles at those same stages include Amh, Csrnp1, Dmrt1, Fdxr, Kdm6b, Pcsk6 and Sox9. Elements with common expression patterns in Chrysemys and Trachemys are candidate key TSD elements for the Emydidae family of turtles to which they belong. A few genes were also differentially transcribed in Apalone late in development, and represent either genes that are thermosensitive in turtles but not involved in TSD, or TSD elements that retained relic thermosensitive expression during GSD evolution. Such genes include Des, Eif4a2 (both upregulated at 31°C during stage 22 in Apalone), Dpt and Twist1 (upregulated at 31°C during stage 19).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Comparative transcriptional patterns in Chrysemys, Apalone and Trachemys for a subset of common elements in the regulatory network of sexual development across embryonic stages at male- (low) and female- (high) producing temperatures for the TSD taxa.

https://doi.org/10.1371/journal.pone.0172044.g006

Clues to TSD from the comparison of turtle versus crocodilian transcriptomes

Surprisingly few differentially-expressed genes were found in common in the categories of interest between Chrysemys (this study) and alligator transcriptomes [24], but the ones detected provide an important insight. Namely, during the TSP of both species Sf1, Amh and Amhr2 were upregulated at MPT, whereas Crabp2 and Hspb1 were upregulated at FPT. Thus, these genes display a sex-specific expression pattern, i.e., upregulation at MPT or FPT across taxonomic groups regardless of relative temperature given that in Chrysemys colder temperatures induce males and warmer values induce females, whereas the opposite is true in alligator [the colder values used in [24] produce females and the warmer temperature produces males]. On the other hand, the heat shock gene Hsph1 was upregulated at MPT in alligator and at FPT in Chrysemys. Hsph1 therefore, displays a temperature-specific pattern of expression (upregulation at warmer temperature in both species) irrespective of the sex-produced. These observations point to Sf1, Amh and Amhr2, Crabp2 and Hspb1 as key common elements in reptilian-TSD for sex-specific development. This finding is consistent with previous research that identified them as critical members of the vertebrate regulatory network for gonadogenesis, and critical during TSD and GSD evolution [e.g., Sf1 [28, 42], Amh and Amhr2 [122–124], Crabp2 [125, 126].

Detection of temporally co-expressed gene clusters

Genes of interest (described in Table 5) clustered during embryogenesis by their co-expression patterns in both turtles. Chrysemys differ more in the number of coexpressed modules (45 modules at 26°C; 10 at 31°C) than Apalone (16 modules at 26°C, 21 at 31°C) (Fig 5), a pattern similar to core eukaryotic genes used as negative control. This suggests that temperature differentially orchestrates gene co-expression in TSD versus GSD, such that Chrysemys’ response is more compartmentalized, and Apalone’s is broader. However, tests with sexed Apalone embryos are needed to rule out any sampling effects that might have contributed to these differences. Some vertebrate sex determination/differentiation genes were clustered in turtles (such as Cbx2 and Dmrt2 in Chrysemys at both temperatures, and Ar/Lhx9 and Insr/Srd5a1 in Apalone), whereas no associations among these genes are known in mammals or birds. Future functional assays on these candidates are warranted. Cluster composition differed by temperature and species, as clusters were enriched in different biological pathways, reflecting temperature effects on gene co-expression and the existence of modules in the urogenital network. Chrysemys male transcriptomes were enriched for pathways regulating transcription, cell proliferation, reproductive development and amino acid phosphorylation (a kinase activity that has been linked to Sry regulation in mice [99, 100]). Immune response functions like lymphocyte and leukocyte activation were female-bias concordant with humans [127]. Cell proliferation, which showed thermosensitive responses here (Table 6, S10 Table), is linked to mammalian sexual development as it is affected by Sry and MAPK signaling [128].

Novel transcripts

A high percentage of novel transcripts in Chrysemys (53%) and Apalone (54%) are currently uncharacterized in SwissProt (which contains manually curated, non-redundant eukaryotic protein sequences) and the ncRNA databases, corroborating genes that were annotated as “predicted” in silico in the Chrysemys genome [51]. Many novel transcripts are male-biased at stages 9, 19 and 22 of Chrysemys (S1 Fig). In conjunction with the greater number of co-expressed clusters discovered at 26°C, this discovery of higher number of novel Chrysemys transcripts at 26°C is indeed curious, and the function of these transcripts merits further investigation.

New insights on known turtle regulators previously studied in GSD turtles

Our RNAseq data shed new light on some important candidate genes previously studied in turtles as follows.

Wt1 regulates the formation of the bipotential gonad, and the maintenance of Sertoli cells and seminiferous tubules in developing testis [129]. Wt1 is upregulated at low temperature before and during the thermosensitive period (TSP) in Chrysemys [42], and in Apalone mutica (GSD) (relic thermal sensitivity in GSD turtles) ([27] and this study) at stage 22 (Fig 5). In contrast, Wt1 expression in mice and chicken gonads is sexually monomorphic through embryogenesis [21, 130]. Finding differential expression from stage 19–22 Apalone gonads is important because previous studies in A. mutica used adrenal-kidney-gonad complexes (AKGs) [27, 28, 31, 34, 35], and expression from the adrenal-kidney can mask gonadal expression, as occurs in Chrysemys and other turtles [42, 131–133].

Sf1 is required for gonadal and adrenal gland formation and steroidogenic activity [134]. Sf1 is directly activated by Wt1, and its expression is thermo-insensitive in GSD turtles ([28] and this study) as in in Chelydra (TSD) [41], but male-biased in Chrysemys ([42] and this study) as in the slider turtle Trachemys scripta (TSD) [135]. This underscores that Sf1 expression is evolutionarily labile, as observed across major vertebrate lineages (reviewed in [42]).

Dax1 is involved in mammalian ovarian and testicular formation [136, 137]. Dax1 is upregulated at low temperature in Chrysemys ([42] and his study), and in Apalone (this study), as observed in A. mutica [31] perhaps driven by the relic thermosensitive expression of its activator (Wt1). In contrast, Dax1 expression is female-biased in birds and monomorphic in several TSD taxa including the green sea turtle Lepidochelys olivacea, Chelydra and Trachemys (reviewed in [42]).

Sox9 is activated by Sry in eutherian mammals, tipping the bipotential gonad towards the male fate. Sox9 is upregulated in Chrysemys ([42] and this study). In contrast, Sox9 in Apalone shifts from upregulation at high-temperature (stage 15) to upregulation at low-temperature (stage 19), perhaps reflecting the evolutionary drift in GSD turtles from its ancestral thermal response. Consistently, Sox9 expression in A. mutica is monomorphic [34].

Aromatase is key in ovarian formation and steroidogenic activity [138, 139]. Aromatase is upregulated at FPT during Chrysemys’ late-TSP ([42] and this study). Aromatase expression is evolutionarily labile across vertebrates [140, 141]. The monomorphic aromatase transcription in the bipotential gonad (stages 9 and 12) in turtles observed here is consistent with its expression in chicken [23].

Dmrt1 regulates sexual development in vertebrates [56] and its molecular evolution is associated with transitions in sex determination in reptiles [126]. Dmrt1 is sex-linked in fish [142] and in birds [143]. Dmrt1 is upregulated in gonads at MPT in Chrysemys (this study; and counter to AKG expression [34]), underscoring Dmrt1’s critical role in testicular development as in Trachemys (TSD) [37]. Dmrt1 in Apalone shifts from low-temperature upregulation (stage 19) to high-temperature upregulation (stage 22) in gonads (this study), whereas its expression is monomorphic in AKG of A. mutica [34], perhaps also reflecting drift of gene expression during GSD evolution as suggested for Sox9 [34]. Dmrt1 expression is male-biased in fish, birds and mammals [23, 42].

4. Are there unique reptilian regulators not present in the gonadal network of other vertebrates?

Discovery of new elements in the vertebrate gonadal network

Our dataset also revealed genes expressed in turtle gonads but unreported in mice, thus possibly unique to reptilian (or turtle) gonadogenesis. Among these are Calr (female-biased at stage 9) and Dcn, a component of the extracellular matrix uncharacterized in the mouse gonad [144], and which showed female-bias in Chrysemys stage 15 onwards. These, as well as the differentially-expressed ncRNAs and unannotated novel transcripts identified here (S2 Table, S1 Fig) in both turtles, represent novel candidates for further study. Some of the differentially-expressed ncRNAs identified here are particularly intriguing, such as 7SK RNA and those classified as putative conserved noncoding region given that these ncRNAs can play key roles in transcriptional regulation [87, 145]. The observed higher expression at FPT of Xist RNA during the TSP of Chrysemys (stage 19), and its absence at other stages and in Apalone is also intriguing. Namely, this is the first report of Xist expression during gonadal development in turtles lacking sex chromosomes [45], a ncRNA region that is critical for dosage compensation via epigenetic inactivation of the X chromosome in human or its upregulation in Drosophila [87, 146].

Finally, genes enriched in hypoxia tolerance and mitochondrial functions that mediate the adaptation to sub-zero temperatures were differentially-expressed in Chrysemys across all stages [147], including translocases that function as chaperones across the mitochondrial membranes, the SLC25 family of mitochondrial transporters [148], Ep300, Casp1 and Thioredoxin family involved in hypoxia signaling [149, 150]. Our data indicate that these genes, which underlie thermal adaptation, are also active from early development.