Transcriptomic responses to environmental temperature by turtles with temperature-dependent and genotypic sex determination assessed by RNAseq inform the genetic architecture of embryonic gonadal development

Vertebrate sexual fate is decided primarily by the individual’s genotype (GSD), by the environmental temperature during development (TSD), or both. Turtles exhibit TSD and GSD, making them ideal to study the evolution of sex determination. Here we analyze temperature-specific gonadal transcriptomes (RNA-sequencing validated by qPCR) of painted turtles (Chrysemys picta TSD) before and during the thermosensitive period, and at equivalent stages in soft-shell turtles (Apalone spinifera—GSD), to test whether TSD’s and GSD’s transcriptional circuitry is identical but deployed differently between mechanisms. Our data show that most elements of the mammalian urogenital network are active during turtle gonadogenesis, but their transcription is generally more thermoresponsive in TSD than GSD, and concordant with their sex-specific function in mammals [e.g., upregulation of Amh, Ar, Esr1, Fog2, Gata4, Igf1r, Insr, and Lhx9 at male-producing temperature, and of β-catenin, Foxl2, Aromatase (Cyp19a1), Fst, Nf-kb, Crabp2 at female-producing temperature in Chrysemys]. Notably, antagonistic elements in gonadogenesis (e.g., β-catenin and Insr) were thermosensitive only in TSD early-embryos. Cirbp showed warm-temperature upregulation in both turtles disputing its purported key TSD role. Genes that may convert thermal inputs into sex-specific development (e.g., signaling and hormonal pathways, RNA-binding and heat-shock) were differentially regulated. Jak-Stat, Nf-κB, retinoic-acid, Wnt, and Mapk-signaling (not Akt and Ras-signaling) potentially mediate TSD thermosensitivity. Numerous species-specific ncRNAs (including Xist) were differentially-expressed, mostly upregulated at colder temperatures, as were unannotated loci that constitute novel TSD candidates. Cirbp showed warm-temperature upregulation in both turtles. Consistent transcription between turtles and alligator revealed putatively-critical reptilian TSD elements for male (Sf1, Amh, Amhr2) and female (Crabp2 and Hspb1) gonadogenesis. In conclusion, while preliminary, our data helps illuminate the regulation and evolution of vertebrate sex determination, and contribute genomic resources to guide further research into this fundamental biological process.


Introduction
Organisms vary wildly in how they determine sex [1]. Vertebrate sex-determining mechanisms range between genotypic sex determination (GSD) and environmental sex determination (ESD) [2,3]. The most common ESD mechanism in vertebrates is temperature-dependent sex determination (TSD). The commitment of the bipotential gonad to differentiate into testes or ovaries is triggered by the genotype in GSD, and by temperatures experienced during the thermosensitive period of embryonic or larval development in TSD [4][5][6]. All studied mammals, birds and amphibians exhibit GSD, while TSD is found in some species of fish and lizards, tuatara, all crocodilians, and most turtles [2,7]. GSD represents systems of extreme developmental canalization while TSD is a textbook example of phenotypic plasticity. Intermediate mechanisms also exist where the genotype is overridden by temperature [2,8], even in the presence of sex chromosomes [9,10]. TSD is ancestral to turtles, reptiles, and perhaps even amniotes [11][12][13][14] such that the TSD of turtles and crocodilians may be considered homologous traits [11,14]. This diversity has defied scientific explanation [1], impeding a full understanding of sex ratio evolution and consequently of population dynamics and diversification both in the past and in the face of contemporary environmental change [15][16][17][18]. This gap is due partly to our incomplete understanding of the molecular basis of GSD and TSD. For instance, the key genetic elements that mediate the effect of environmental temperature in TSD systems remain elusive. Unlike GSD models such as mammals [19][20][21][22] and chicken [23] whose gonadal developmental pathways are well understood (albeit not fully), our knowledge for GSD reptiles and TSD species is incipient, and prevents deciphering the genetic architecture of vertebrate sex determination and its evolution.
The first reptilian gonadal transcriptomes were recently characterized in alligator (TSD) [24], and slider turtles (TSD) [25], whereas previous studies on GSD and TSD turtles had used exclusively quantitative PCR and in situ hybridization targeting a number of genes underlying sexual development (full gene names are found in Table 1 Table). While these earlier turtle studies provided only fragmentary information [30,33,41,42], their use of sensitive qPCR permitted the detection of differential gene expression that can pass undetected in transcriptomic analyses [43,44], and these efforts uncovered substantial evolution of transcriptional patterns for some elements across vertebrates. However, the extent of evolution of the composition and expression of this gene regulatory network remains unclear. Thus, deciphering the composition, environmental sensitivity, and evolution of the gonadal gene network in additional TSD turtles and how they compare to GSD turtles is overdue.
Here, we use a comparative approach to test for transcriptional responses to environmental temperature (or lack thereof) at several stages of embryonic development in two turtle species, the painted turtle Chrysemys picta (TSD) and the soft-shell turtle Apalone spinifera (GSD), hereafter denoted as Chrysemys and Apalone respectively. Chrysemys is a TSD turtle lacking sex chromosomes [45] whose thermal ecology has been studied extensively [46][47][48][49], and constitutes an emerging model for ecology, evolution, and human health [50] with increasing

Fhl2
Four And A Half LIM Domains 2

Fog2
Friend of Gata Protein 2

Gata2
Gata binding protein 2 genomic resources available [51][52][53]. Apalone is a GSD turtle with ZZ/ZW sex chromosomes [54] whose sex ratios are unaffected by the incubation temperature [55], and thus, it serves as a negative control for TSD responses.
Our RNAseq approach provides the first glimpse of the full transcriptional network in closely related reptiles with contrasting sex-determining mechanisms (GSD and TSD), and complements recent contrasts between turtles and model mammalian developing gonads [25], as well as enabling a comparison with alligator (another TSD reptile) [24]. With these data we carry out an initial test of several hypotheses to help elucidate whether the transcriptional building blocks in TSD and GSD systems are identical, and whether differential deployment of common elements distinguishes these systems. First, we examine whether the vertebrate sex determination/differentiation gene regulatory network known from model mammals and birds [3,42,56] is present and active in painted turtles. Second, we test whether the conversion of thermal inputs into TSD sex-specific development in painted turtles might involve the epigenetic machinery [57], hormonal pathways [58], or general sensing responses [59], by examining the transcriptional response to incubation temperature of heat-shock genes, transient receptor potential genes, germ-line and histone-related genes, as well as androgen-and estrogen related genes. Third, we test for differences in the transcriptional response to temperature between TSD and GSD turtles to uncover changes that may set these systems apart. Additionally, we identify novel candidate genes in painted turtles undescribed in mice and compare this information to that of other TSD reptiles [24,25] to test whether they represent potentially unique reptilian regulators. Importantly, our transcriptomic time series of developing gonads at high and low temperatures in TSD and GSD turtles across developmental stages before, during and after the activation of the thermosensitive period (TSP) constitutes a critically needed resource to facilitate more extensive research to illuminate the proximate ecological regulation and evolution of vertebrate sex determination under various thermal regimes [16]. We note that for each species a single transcriptome was obtained from pooled embryos at each developmental stage at each temperature, such that our results are preliminary, and our conclusions represent critical working hypotheses to foster further research in this fascinating field.

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-kidneygonadal (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 NanoDrop1 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][76][77]. Briefly, Fisher's exact test is based on a 2x2 contingency table (Table 2) 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 11 + n 12 is the marginal total for focal gene X; n 21 + n 22 is the marginal total for the remaining genes in the transcriptome; n 11 + n 21 is the marginal total for the 26˚C treatment; n 21 + n 22 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 0 : θ = 1, where θ = π 11 π 22/ π 12 π 21 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 differentiallyexpressed 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-kidneygonad 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; Dmrt1exon 1 is mis-assembled/mis-annotated; Sox9-exon 1 is mis-assembled).
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.

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   Table 5. Differential expression was assessed between temperatures at each developmental stages for each species separately. doi:10.1371/journal.pone.0172044.g001 Transcriptomic response to incubation temperature in developing gonads of TSD and GSD turtle embryos 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   Table).
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.

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.

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 RNAsplicing) 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).

Discussion
The evolution of sex determination remains an evolutionary enigma, yet this developmental process is critical for the production of sex ratios and consequently for the dynamics and evolutionary potential of populations [92]. Turtles are an ideal study system since TSD and GSD co-occur in this group, yet relatively little is known about the gene network controlling turtle urogenital development. Here we present an initial characterization of the full composition of this network using RNAseq and its transcriptional response to incubation temperature in 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]. doi:10.1371/journal.pone.0172044.g003 developing gonads of TSD and GSD turtles. This effort sheds light on the genome-wide architecture of vertebrate gonadogenesis and the evolution of its environmental sensitivity in TSD and GSD systems [45][46][47][48][49][50][51][52][53][54]. While the expression of some candidate genes has been profiled in a few TSD reptiles (S1 Table) and two transcriptome analyses were recently reported in TSD turtle and alligator [24,25], Apalone remains to our knowledge the only turtle GSD genus [55] whose embryonic urogenital transcription has been studied (this study and [27, 28, 31, 34, 35]. We note that because molecular sexing was unavailable for Apalone when data were collected for this study, only temperature effects (and no sex effects) could be analyzed here. Future studies will be able to leverage novel molecular sexing techniques [93] for a more comprehensive test of temperature, sex, and their interaction on transcription in Apalone. Overall, our data provides circumstantial evidence that the transcriptional circuitry underlying gonadogenesis in TSD and GSD turtles is broadly the same, and that differences between these mechanisms are largely due to the differential deployment of these common elements as detailed below. We identified novel candidate genes whose early differential expression suggest that they may contribute to transmitting the temperature signal to the developmental pathway, potentially helping determine the sexual fate in TSD turtles, as well as candidate genes undescribed in the gonadal regulatory network of mice and chicken. We first discuss general characteristics of the transcriptomes generated here and then address a series of questions related to our central hypotheses.

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)

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][95][96][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.

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 androgenand 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 histonerelated 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 temperaturespecific 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][20][21][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.

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][17][18][19][20][21][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). 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][123][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][132][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

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 differentiallyexpressed 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.

Conclusion
Ours is the first comparative transcriptomic analysis of TSD vertebrates and GSD turtles. Our analyses uncovered numerous homologs of mammalian urogonadal genes that were unidentified in turtles, alongside many turtle-specific novel transcripts. The strengths of the transcriptomic time series through embryogenesis in two species permitted an initial and simultaneous characterization of the genome-wide network underlying gonadogenesis and its response to environmental temperature in TSD and GSD systems. Because a single pool of RNA was available per temperature-by-stage, our inferences should be taken as a strongly suggestive initial glimpse of the genome-wide composition and regulation of the gene network underlying sexual development in TSD and GSD turtles. These results, accompanied by the discovery of previously unknown gonadal regulators in vertebrates that are active well before the onset of the thermosensitive period highlight the value of targeted investigations of early orchestrators of embryogenesis. Further, the thermosensitive response of key elements of multiple signaling pathways potentially governing sex determination observed here, underscore that differences between TSD and GSD in turtles are less likely due to unique elements in this network (although some candidates were identified) but instead, perhaps largely due to the differential deployment of common elements and modules. Our work thus contributes to deciphering the evolutionary puzzle of vertebrate sex determination, and provides significant genomic resources and working hypotheses to guide further investigations in this active field of research.  Table. Comparisons of genes of interest (described in Table 5 Table. Genes whose embryonic gonadal expression was explored in this study (symbols and names).