Gene expression of male pathway genes sox9 and amh during early sex differentiation in a reptile departs from the classical amniote model
BMC Genomics volume 24, Article number: 243 (2023)
Sex determination is the process whereby the bipotential embryonic gonads become committed to differentiate into testes or ovaries. In genetic sex determination (GSD), the sex determining trigger is encoded by a gene on the sex chromosomes, which activates a network of downstream genes; in mammals these include SOX9, AMH and DMRT1 in the male pathway, and FOXL2 in the female pathway. Although mammalian and avian GSD systems have been well studied, few data are available for reptilian GSD systems.
We conducted an unbiased transcriptome-wide analysis of gonad development throughout differentiation in central bearded dragon (Pogona vitticeps) embryos with GSD. We found that sex differentiation of transcriptomic profiles occurs at a very early stage, before the gonad consolidates as a body distinct from the gonad-kidney complex. The male pathway genes dmrt1 and amh and the female pathway gene foxl2 play a key role in early sex differentiation in P. vitticeps, but the central player of the mammalian male trajectory, sox9, is not differentially expressed in P. vitticeps at the bipotential stage. The most striking difference from GSD systems of other amniotes is the high expression of the male pathway genes amh and sox9 in female gonads during development. We propose that a default male trajectory progresses if not repressed by a W-linked dominant gene that tips the balance of gene expression towards the female trajectory.
Further, weighted gene expression correlation network analysis revealed novel candidates for male and female sex differentiation.
Our data reveal that interpretation of putative mechanisms of GSD in reptiles cannot solely depend on lessons drawn from mammals.
Gonads are unusual in that two different organs, testis or ovary, develop from the same bipotential primordium, the genital ridge. In species with genetic sex determination (GSD), the fate of the bipotential gonad is determined at the point of conception by the complement of chromosomes inherited from the parents. Mammals usually have GSD with male heterogamety (XX females/XY males) whereby the presence or absence of the Y chromosome, and the gene SRY it harbours, determines the developmental trajectory leading to the sexual phenotype. SRY initiates the male regulatory pathway which leads to suppression of the female pathway. Without SRY the female regulatory pathway is pursued, which leads to suppression of the male pathway . Birds also have GSD but with female heterogamety (ZW females/ZZ males). Here, dosage of DMRT1, a gene encoded and functional on the Z but not the W chromosome, determines whether the female or male regulatory pathway is taken [2, 3]. Only a two-fold difference in expression of DMRT1 (expressed from two Z chromosomes in males versus only one Z in females) is sufficient to determine sexual fate. The two-fold expression difference occurs at a blastoderm stage before the mesonephros or genital ridges form, but in the gonadal lineage, expression of dmrt1 is highly upregulated maintaining the two-fold dosage difference that determines sexual fate of the gonads . In both birds and mammals, the sex determining gene (DMRT1 or SRY, respectively) initiates the run-away upregulation of SOX9, a key transcriptional activator of male sex genes, a process required for testes differentiation . With respect to sex determination, fishes are an interesting and the most diverse group, as all so far discovered vertebrate sex determination systems can be found among them .
Although the mammalian and avian GSD systems have been well studied, in reptiles no female or male GSD pathway has been investigated. Part of the reason is that, unlike mammals and birds, reptiles display a variety of sex determination systems, including genetic sex determination (GSD) and temperature sex determination (TSD) or a combination of the two. For instance, the central bearded dragon, Pogona vitticeps, has a ZZ male: ZW female system, but shows male to female sex reversal at elevated temperatures. Some reptiles with GSD have highly differentiated sex chromosomes (either male or female heterogamety), whereas others have poorly differentiated or homomorphic sex chromosomes . Studies of the GSD turtles Apalone spinifera and Apalone mutica have been published but no female or male sex differentiation network has been described [8,9,10,11]. An efficient method to determine the sex chromosome complement in embryos is crucial for studying male and female GSD pathways, but was lacking for A. spinifera and A. mutica in the past . Such sex chromosome PCR markers have been long established in Pogona vitticeps , the first reptile model species to have emerged for focused studies of male and female GSD pathways.
Most work on reptile sex determination pathways has concentrated on TSD species. In TSD the egg incubation temperature during the temperature sensitive window of development biases the percentage of male or female offspring . In some species, the higher temperature leads to females, in others to males . The environmental temperature signal is captured by ancient and ubiquitous signalling pathways in the cell that influence the epigenetic release or suppression of key sex genes [16, 17] to direct developmental programming of sex in the embryo [18, 19]. Downstream genes in the sex determining network of TSD species are similar to those acting in mammalian or avian GSD. For example, dmrt1, sox9 and amh are involved in the male trajectory, foxl2 and aromatase (cyp19a1) in the female trajectory [20,21,22].
We therefore set out to discover how gene expression profiles differ between male and female developmental trajectories under the control of the sex chromosomes in Pogona vitticeps and determine when in embryonic development sex differences first emerge. We compared expression profiles in the gonads of embryos of Pogona vitticeps incubated at moderate temperature (28 °C) at which no sex reversal occurs [13, 26, 27] (Fig. 1A). We demonstrate that sex specific differences have already emerged by the earliest stage at which we can dissect a consolidated gonad (stage 6 [28, 29]). Thus, sex determination under GSD occurs at a very early stage of gonadal differentiation. Patterns of expression of key genes in the vertebrate male sex determining network, sox9 and amh, differ strikingly from those established in mammalian models. We found early engagement of amh in the developing testes. However, early sex-nonspecific expression of sox9 leads us to propose that SOX9 is not a central hub for early reptile sex differentiation and there may be no positive feedback loop driving up the expression of SOX9 specifically in the male developmental trajectory of reptiles as it is of mammals and birds.
For the purpose of this paper, sex determination refers to the processes by which sexual fate is decided and sex differentiation refers to the regulatory processes that lead to development of either a testis or ovary. Sex determination in GSD species is most commonly associated with differential expression of genes encoded on sex chromosomes.
Male and female gonad transcriptome profiles diverge at an early developmental stage
We compared gene expression of female and male gonads at three developmental stages, when gonads begin to differentiate (stage 6/7), are at an early stage of differentiation (stage 12) and are differentiated (stage 16), as well as in adults (Fig. 1A) [28, 29].
Female and male gonadal transcriptomes were distinguishable at the earliest stage at which we can dissect a consolidated gonad,stage 6/7 (Figure S1A); 222 genes were expressed differentially in male and female embryos at stage 6/7 (differentially expressed genes; DEGs; adjusted p-value < 0.05), with 130 of those genes showing a male bias and 92 showing a female bias (Fig. 1B, D; Table S2). As expected, over the course of ovary and testis development the number of DEGs increased (to ca 1250 at stage 12 and 16). Roughly the same proportion of genes was more highly expressed in ZZ or ZW developing gonads (Fig. 1B, D). As development progressed, the transcriptomic profiles of ZZ males and ZW females diverged further, as reflected by the increase in number of DEGs to 8,425 between adult ovaries and testes (Fig. 1D and S1B).
During embryonic development, genes typically had either a male or female bias; only 3 exchanged roles (chrna1, Cholinergic Receptor Nicotinic Alpha 1 Subunit; noxo1, NADPH oxidase organizer 1; and ENSPVIG00000017292, keratin type II cytoskeletal 5-like). A cluster heatmap of the 300 most relevant DEGs (lowest adjusted p-value) showed that one set of genes was upregulated in females and another set of genes was upregulated in males, each with a different degree of repression in the other sex (Fig. 1E). However, a few genes had the opposite bias during development and in adulthood (Figure S1E), presumably fulfilling different roles in the embryo and the adult.
Roughly half of the genes (52%) maintained their differential expression levels during embryonic development, once their expression had diverged (Fig. 1C and S1C, D and Table S2 columns J, K). However, 106 DEGs were specific for stage 6/7 and 554 DEGs were specific for stage 12 (Fig. 1C). Of the 106 DEGs that were dimorphic only at stage 6/7, 70 were more highly expressed in females and 36 were more highly expressed in males. In order to establish equal expression between males and females by stage 12, gene expression for most of those genes changed in female embryos adjusting to the expression in male embryos (Figure S2A). Those 106 stage 6/7 specific DEGs were not typically implicated with sex differentiation. The male biased group showed no gene ontology enrichment. However, the female biased group was enriched for gene ontology terms related to muscle development/contraction (Figure S2B). The 21 driver genes of that enrichment were much more highly expressed in female than in male gonads (4 to 64-fold), implying biological relevance. Expression of all but one of these genes was progressively reduced during female development until, by stage 12, they matched the low expression level of male gonads (Figure S2C-F). None of the 21 muscle related genes have been previously implicated in gonadal development or sex differentiation. It is difficult to see why muscle contraction should be implicated in ovary development, and the role of several of these genes in calcium regulation may be more pertinent, as calcium signalling is involved in sex reversal in P. vitticeps [18, 19]. Amongst the group are genes involved in Ca2+ signalling such as myoz2 (Myozenin 2) which tethers calcineurin , adrb2 (β-2-Adrenergic Receptor) which is associated with class C L-type calcium channel , subunits of troponin (tnni1, tnni2, tnnc2, tnnt2, tnnt3) which regulate muscle contraction via Ca2+ binding , and subunits of the nicotinic acetylcholine receptors (chrnd and chrng) which are non-selective cation channels . It will be interesting to investigate if this group of genes is transiently expressed specifically in the female trajectory or if those genes are expressed in the precursor cell line, the intermediate mesoderm which develops into the gonads and kidneys. The switch in transcriptome reprogramming (i.e., the downregulation of redundant genes) might simply happen faster in the male trajectory.
Nr5a1 is involved in gonad development but probably not in sex differentiation
Two genes that are highly expressed in undifferentiated gonads in both sexes of mammals and chicken are nr5a1 (Nuclear Receptor Subfamily 5 Group A Member 1) and gata4 (GATA Binding Protein 4) [2, 34, 35]. These marker genes were also highly expressed in P. vitticeps in both sexes from the earliest stage examined (Fig. 2A, B, right panels). Both genes continued to be highly, but not differentially, expressed, in both sexes through to stage 16. Adult gonads showed much lower expression of nr5a1 and gata4, suggesting their main role is during embryonic development (Fig. 2A, B). Their similar expression patterns suggest a common role for nr5a1 and gata4 in the early formation of the genital ridge in amniotes.
Nr5a1, encoding the Steroidogenic Factor 1 (SF1), resides on the sex chromosomes of P. vitticeps . Coupled with its role in sex determination in mammals, this location makes it a candidate for the master sex determining gene . However, nr5a1 showed no sex bias in expression levels that might be expected for a gene with a role in sex determination or differentiation.
Of the 302 protein coding genes (including nr5a1) and long non-coding RNA genes borne on one of the four sex chromosome scaffolds of P. vitticeps [24, 36], only one gene showed sex biased expression at stage 6/7, displaying a ca 2.6-fold higher expression in females (adj. p-value = 0.044). This was ak8 (Adenylate Kinase 8). Adenylate kinases are small enzymes with a crucial role in energetic metabolism. They are nucleoside monophosphate kinases that catalyse the reversible transfer of the terminal phosphate group between nucleoside triphosphates and monophosphates (Figure S3A) . Isoenzyme 8 is cytosolic and least studied of all isoforms . As of now, it has not been involved in sex determination or differentiation. Another 29 sex chromosome genes (ca 10%) were differentially expressed between the sexes at later stages 12 and/or 16 (Figure S3A).
Chromatin remodelling genes jarid2 and kdm6b are not differentially expressed in GSD
Histone modifiers jarid2 (Jumonji And AT-Rich Interaction Domain Containing 2) and kdm6b (Lysine Demethylase 6B) play a role in high temperature induced sex reversal in P. vitticeps [19, 38] and in other evolutionary disparate TSD species [16, 17, 20, 39,40,41]. Epigenetic regulation is also implicated in GSD sytems, for example in mice, where a Jumonji family member, Jmjd1a (Lysine Demethylase 3A), directly regulates expression of the sex determining gene Sry .
We found that during GSD in P. vitticeps both genes, jarid2 and kdm6b, were highly expressed but showed no sex bias (Figure S4 A, B). Two transcript isoforms of jarid2 and kdm6b are expressed in P. vitticeps ZZ and ZW embryonic gonads, but no significant difference was observed in their ratio between ZZ and ZW gonads (Whiteley et al. 2022, Science Advances).
Of 154 chromatin modifying genes (PANTHER protein class PC00077), only one was differentially expressed in P. vitticeps at embryonic stage 6/7, though at low abundance. Samd7 (Sterile Alpha Motif Domain-Containing Protein 7) is a retina specific component of the chromatin modifier complex PRC1  (Table S6). Only a few other chromatin modifying genes were differentially expressed between males and females at later embryonic stages, whereas roughly half of the chromatin modifying genes were differentially expressed between adult ovaries and testes (Table S6). This suggests that epigenetic regulation plays a major role in maintenance of differentiated adult gonads rather than in establishing differentiation during embryonic development.
Known sex differentiation genes have unique expression profiles in P. vitticeps
Among transcription factors that were highly expressed in both sexes at stage 6/7 (Table S3) was sox9 (SRY-Box Transcription Factor 9), a gene with a critical role in vertebrate sex determination. Sox9 is autosomal in P. vitticeps [24, 36], implying that it is not the master sex determining gene in this species. We found that sox9 expression was high in both sexes at stage 6/7. High expression is maintained in males through stage 12 and 16 but decreased with time in females (Fig. 2C). This pattern of expression of sox9 in P. vitticeps is quite different from that in mammalian and avian systems, in which an initially low expression of SOX9 is dramatically upregulated in males in response to differential expression of the sex determining gene [2, 44] (Fig. 3A).
FGF9 (Fibroblast Growth Factor 9) and its receptor FGFR2, are highly expressed in male gonads and positively regulate SOX9 expression during testis development in mammals [45, 46]. However, in P. vitticeps the expression of these genes showed no male bias at the beginning of differentiation, stage 6/7, and at stages 12 and 16 fgf9 expression was higher in females than in males (Fig. 2D, S4C). Expression profiles of sox9 and fgf9 showed no correlation. In males sox9 expression from stage 6/7 to stage 16 was constant while fgf9 expression decreased. In females fgf9 expression was constant while sox9 expression decreased (Fig. 2C and S4A). Thus the positive feedback loop between sox9 and fgf9 in mammals seems not to operate in P. vitticeps, at least after stage 6/7. However, it remains possible that very early in development (prior to stage 6/7) sox9 upregulation is supported by fgf9 in both, male and female P. vitticeps, then later a dominant gene may suppress sox9 expression in females despite constant levels of fgf9 (Fig. 3A).
AMH (Anti-Müllerian Hormone) is differentially expressed early in mammalian sex determination, being upregulated by SOX9 in male gonads and absent in female . We found that amh (ENSPVIG00000005953) is highly expressed in gonads at stage 6/7 in P. vitticeps females, and even more highly (about fourfold) expressed in males (Figs. 2E and 3B). Expression of amh remains constant in females, but in males increases to 16-fold over females during embryonic development. In adults of both sexes it is reduced, but still sexually dimorphic (Fig. 2E).
Dmrt1 (Doublesex And Mab-3 Related Transcription Factor 1) is early acting in birds , as befits its role as the bird sex determining gene, but it has a later action in testis differentiation in mammals, especially mice . In P. vitticeps, dmrt1 is differentially expressed early in gonad development, showing a ca 19-fold higher expression in ZZ than in ZW gonads at stage 6/7, one of the strongest male biases early in development (Fig. 1B). In males, expression of dmrt1 stayed high during the remainder of testis development through to stage 16. In females expression of dmrt1 increased over time, from ca 19-fold lower expression than the male gonad to only ca fourfold lower expression by stage 16 and in adult ovaries (Figs. 1B and 2F). This expression profile strongly indicates an important role for dmrt1 in early gonadal sex differentiation in P. vitticeps but it cannot be the master sex regulator due to its location on autosomal chromosome 2 and not the Z/W sex chromosomes .
In mammals, WNT signalling is active in the undifferentiated gonad, mediated by positive regulators WNT4 (Wnt Family Member 4) and RSPO1 (R-Spondin 1); after SRY activation in the male gonad, WNT4 and RSPO1 expression is downregulated to repress the female pathway [49,50,51]. In P. vitticeps, wnt4 and rspo1 were well expressed throughout embryonic development in both sexes; neither gene was downregulated in males (Fig. 4A, B). In fact, no positive regulators of WNT signalling (WNT ligands and frizzled receptors, including wnt4 and rspo1) were differentially expressed at stage 6/7 (Table S4). Instead, five genes known to antagonize WNT signalling (frzb, trabd2a also known as tiki1, shisa8, dkk3 and wif1)  all showed a male bias at stage 6/7 of gonad development (Figure S4D-H, Table S4). This suggests that control of WNT signalling in male P. vitticeps early in gonad differentiation is achieved by expression of WNT inhibitors rather than by suppression of WNT activators wnt4 and rspo1. After the onset of differentiation, the expression of WNT pathway genes was more diverse. Several WNT ligands (wnt2b, wnt4, wnt10a and wnt11) displayed a female bias, while wnt6 expression was male biased. Likewise, several WNT antagonists were male or female biased at stages 12 and/or 16 (Table S4).
The transcription factor FOXL2 (Forkhead Box L2) is expressed early and is female specific in mammalian gonads [53, 54]. In all female vertebrates FOXL2 is likely to be a key regulator of expression of CYP19A1 (aromatase) , the enzyme controlling estrogen production . Estrogen has a highly conserved role in nonmammalian vertebrates, but its function in the development of the mammalian ovary remains less clear . Female markers cyp19a1 and foxl2 were more highly expressed in P. vitticeps females compared to males throughout all embryonic stages (Fig. 4C, D). Both genes were barely expressed in males. While foxl2 expression in females had reached its height by stage 6/7 and was constant throughout development, expression of cyp19a1 rose by roughly eightfold from stage 6/7 to 12, in line with the established view that foxl2 regulates the expression of aromatase .
Novel candidates for sex differentiation
Many other genes were dimorphically expressed throughout all embryonic stages; 17 with a female bias and 78 with a male bias (Table S5). Gene ontology terms enriched in both groups relate to sex development. Female specific terms included tube/vessel regulation and hormone secretion/metabolism. Male specific terms included regulation of meiosis, signalling, development and cell death (Figure S5A, B).
The five most highly ranked genes in the male biased group with at least 16-fold sex difference were amh, col9a1 (Collagen Type IX Alpha 1 Chain, ENSPVIG00000019420), nkx3-1 (NK3 Homeobox 1), fmod (Fibromodulin) and spink2 (Serine Peptidase Inhibitor Kazal Type 2, ENSPVIG00000008173) (Table S5, Figure S5C-F). Of these genes only amh has a known role in embryonic gonad development . NKX3-1 is a transcription factor expressed in the mammalian prostate and throughout its development, and is a well-established prostate tumour suppressor [58, 59]. A recent study demonstrated NKX3-1 expression in human Sertoli cells by immunohistology , although an earlier study detected only trace amounts of NKX3-1 transcripts in human testes . However, Nkx3-1−/− knock-out rats displayed non-aberrant testis . SPINK2 is a protease inhibitor that was shown to have a function in spermatogenesis but has not yet been reported as expressed during embryonic gonad development . Based on their relevant expression pattern for male gonad differentiation in P. vitticeps, both genes are prime candidates for driving male gonad differentiation, in conjunction with established factors, such as dmrt1 or amh. COL9A1 and FMOD are two extracellular matrix proteins. Col9a1 was reported to be expressed with a male bias in the bipotential gonad of an amphibian, Xenopus laevis . In mouse and chicken, the alpha 3 chain, col9a3, has been associated with male gonadal development, but it is unclear if COL9A3 is functioning alone or in a trimer together with the alpha chains 1 and 2 forming collagen type IX . Likewise, the function of collagen type IX in male gonadogenesis is still unknown. The role of cell adhesion or the extracellular matrix (ECM) in establishing testicular or ovarian structure has not been well studied . Our results from P. vitticeps suggest a likely role of ECM genes in early gonad differentiation and provide candidate genes for further studies.
The highest ranked gene in the female biased group was aromatase, cyp19a1 (Table S5). To find more genes relevant for the female trajectory we performed a weighted correlation network analysis that yielded 25 modules (Figure S6A, Table S7). Two modules (module 8/pink and module 14/cyan) strongly correlated with female development (p-value < 0.01, correlation value > 0.7) (Figure S6B). Together they comprised 633 genes, amongst them aromatase cyp19a1, with the highest module-membership in module 14. Several signalling pathways were enriched in this gene set, MAPK, RAS, RAP1, PI3K-Akt, VEGF and calcium signalling (in ascending order of their false discovery rate) (Figure S6C). MAPK and PI3K-Akt pathways have been implicated in female gametogenesis in mammals  and it will be exciting to explore their influence in embryonic female gonad differentiation. Calcium signalling was suggested to be modulated by gonadal hormones  and in P. vitticeps, calcium signalling is thought to be important for the capture of the temperature cue and its transduction to cellular signalling pathways in male to female sex reversal [18, 19]. Our findings suggest that calcium signalling also plays a role in genetically determined female gonad development. VEGF and RAP1 signalling have not been implicated in ovary related functions, though RAP1 plays a role in mouse spermiogenesis . RAP1 is a key regulator in cell junction formation  and VEGF signalling promotes endothelial cell growth, migration and survival, and plays a key role in angiogenesis . The most significant enrichment was in genes involved in the pathway for axon guidance, which presumably have a different function in female gonad development, so might be moonlighting proteins .
Visualisation of the weighted gene correlation network helped to elucidate hub genes in female and male sex differentiation (Figure S7-9, Table S7). Two genes stood out as hub genes with strong positive membership in module 14/cyan (Figure S7A, Table S7), the module with the most significant and highest correlation with ZW genotype. One of them was aromatase cyp19a1. The other, transcription factor lef1 (Lymphoid Enhancer Binding Factor 1), was shown to interact with ß-catenin to regulate target genes of WNT signalling . It has not been strongly associated with gonadal sex differentiation so far. Hints come from two studies showing female-specific regulation of LEF1 splicing in mouse gonadal transcriptomes , and downregulation of lef1 in an ovary-to-testis gonadal transformation process in zebrafish . Being a hub gene in female sex differentiation in P. vitticeps, lef1 might function to orchestrate WNT signalling in female development. The most prominent hub gene in the second ZW module, 8/pink, was hsd17b3 (Hydroxysteroid 17-Beta Dehydrogenase 3), a gene with a strong negative module-membership (Figure S8, Table S7). It is barely expressed in female gonads at any stage but highly expressed in male gonads. Hsd17b3 has a well-known role in male sex differentiation as it catalyses the final step in the testosterone biosynthesis pathway . Suppression of hsd17b3 is important for female development . We are providing with the unsigned correlation network a wealthy source for further investigation to better understand how hsd17b3 is suppressed in female gonads. For example, two transcription factors that are connected with hsd17b3 through a negative correlation were tshz3 (Teashirt Zinc Finger Homeobox 3) and lmx1a-like (LIM Homeobox Transcription Factor 1 Alpha—Like, ENSPVIG00000024375), both, themselves prominent hub genes of module 8/pink with a positive module membership (Figure S8, Table S7). They were more highly expressed in female than male gonads throughout all stages (Table S2). Tshz3 has been recently implicated into Müllerian duct formation in chicken  and lmx1a into ovary morphogenesis in the fruit fly .
Module 20/royalblue with significant correlation with ZW genotype (p-value = 0.006, correlation value = 0.45) comprised the earlier identified group of genes related to muscle development/contraction that was much more highly expressed in female than in male gonads at stage 6/7 and reduced expression until they matched the low expression level of male gonads by stage 12 (Figure S2B-F), This module formed a very concentric and balanced network (Figure S7B).
One module, 7/black, strongly correlated with male development (p-value < 0.01, correlation value > 0.7) (Figure S6B, S9, Table S7). The three most prominent hub genes, all with positive membership, were met (MET Proto-Oncogene, Receptor Tyrosine Kinase) and tiam2 (TIAM Rac1 Associated GEF 2), which have not been associated with sex differentiation, and plcb1 (Phospholipase C Beta 1) which was implicated recently in development of genitalia in male geese . The next three most prominent hub genes were the transcription factor nkx3-1 and the two extracellular matrix genes fmod and col9a1, three of the earlier mentioned genes with strong male bias in all three stages. They all constitute interesting and promising candidates for further exploration.
We undertook the first transcriptome wide analysis of embryonic gonad differentiation in a reptile with genetic sex determination. We found that gonad transcriptome profiles of ZZ male and ZW female embryos diverged early, being already differentiated at stage 6/7 [28, 29]. This suggests that determination of the male or female pathway in P. vitticeps happens before stage 6/7 and thus prior to the consolidation of the gonad into a distinct organ, and possibly before eggs are laid by the mother.
The master sex gene that initiates gonad differentiation is unknown in P. vitticeps and we hoped that it might be revealed by the differential expression of a gene on the sex chromosomes at an early stage. Based on our current understanding of other GSD systems, there are two possible mechanisms by which sex could be genetically determined in P. vitticeps: a female-dominance system similar to the male-dominance system in mammals, and a dosage dependent system similar to that in birds. These might be expected to yield different expression patterns. The mammalian male-dominant sex determiner, the Y-borne gene SRY, is expressed highly for only one day in mice (though more broadly in other mammals), and this is enough to initiate the male pathway [81, 82]. A dosage dependent system, such as the two-fold higher expression of DMRT1 in male chickens , probably requires expression for longer during bipotential gonad development.
The strongest candidate sex determining gene identified so far in P. vitticeps is nr5a1, which has alleles on the Z and W sex chromosomes. However, we observed no differential expression at stage 6/7 of nr5a1, or other sex relevant genes on the sex chromosomes. It remains possible that our sampling missed a short burst of activity of a male- or female-dominant sex determiner before stage 6/7, or that subtle dosage shifts of a dosage-controlled gene were beyond the resolution of this study. Recently, we have described sex-specific transcript isoforms of nr5a1 in gonads of adult P. vitticeps, which are likely to result in NR5A1 proteins with sex-specific function . The technique of high-throughput short read sequencing used in this study is suboptimal to detect those isoforms as they differed in a GC rich and repetitive region. We cannot exclude the possibility that sex-specific isoforms or sex-specific post-transcriptional, translational or post-translational regulation of nr5a1 or another sex chromosome-borne gene may influence sex determination. Genome-scale studies of the translational state to better understand post-transcriptional regulation of sex differentiation genes are needed.
There is also the possibility that other Z or W-borne genes will be discovered that have been missed from the P. vitticeps draft genome assembly [24, 36]. An improved chromosome scale genome assembly, ideally with phased sex chromosomes, may be required to identify the master sex gene in P. vitticeps.
SOX9 is a key gene in mammalian and avian testis determination [5, 84, 85]. It is directly upregulated in males by the sex determining factor, SRY in mammals and DMRT1 in chicken . SOX9 is not upregulated in female mammals. As well as directly initiating AMH expression to effect Müllerian duct regression, SOX9 activates testis-differentiation genes and inhibits ß-catenin in the ovary differentiating pathway . The expression profile of sox9 in P. vitticeps suggests that SOX9 does not act as a central key player in early testis determination in this species, implying that it plays a different role than in mammals. We observed that sox9 was highly expressed in female as well as male embryos at an early stage, then was later, after the onset of gonad differentiation, downregulated in females. Similar sox9 expression was reported in embryonic gonads of the TSD turtle T. scripta [22, 86, 87]. Hence, in the absence of the necessary correlative and functional work on reptiles, interpretation of putative mechanisms of sex determination in reptiles cannot depend on lessons drawn from mammals.
Gonads of T. scripta initially develop morphologically like a testis, with primitive SOX9-positive sex cords that later degenerate at a female producing temperature when sox9 expression is downregulated [22, 86]. We observed no similar development of cord-like structures in female P. vitticeps gonads . The high initial sox9 expression in P. vitticeps female embryos evidently does not masculinize female P. vitticeps embryos, suggesting either the action of a female specific inhibitor or that, in contrast to mammals, SOX9 in is not important in the regulatory pathways that dictate sexual fate in P. vitticeps or reptiles generally.
The function of a female-specific SOX9 inhibitor might be fulfilled by aromatase, an enzyme responsible for a key step in the biosynthesis of estrogens , which is expressed early and strongly in female P. vitticeps embryos. In a marsupial mammal and a human-derived testis cell line, estrogen can regulate gonadal development through the nucleocytoplasmic shuttling of SOX9. In particular, it can supress nuclear translocation of SOX9, so that SOX9 cannot fulfill its function as transcription activator of male pathway genes [89, 90]. A similar mechanism might prevent SOX9 in female P. vitticeps embryos from activating male pathway genes in the nucleus until after sox9 expression decreases later in female embryonic development. Next to estrogen, also the signalling molecule prostaglandin D2 was shown to induce nuclear import of SOX9 in mammals which induces male gonad development [91,92,93]. We cannot draw any conclusion for P. vitticeps, as the gene for the prostaglandin synthase (L-Pgds or H-Pgds) has not yet been annotated in its genome. It will be exciting to explore the cellular location of SOX9 protein in embryonic gonads of P. vitticeps.
DMRT1 and SOX9 have been shown to functionally interact in mammals . The lower expression level of dmrt1 in female gonads of P. vitticeps at the early stage of development might prevent SOX9 from being fully functional.
The regulation of expression of SOX9 transcripts differs in mammals and other vertebrates. In mammals, SOX9 is activated by the Y chromosome borne SRY gene in male (but not female) embryos as an early event in testis determination. Then NR5A1, FGF9 and SOX9 itself establish and maintain high SOX9 expression in male gonads. In chicken, DMRT1 is thought to initiate SOX9 expression. A direct role for DMRT1 seems unlikely in P. vitticeps because dmrt1 and sox9 expression profiles were inversely correlated in females. Although, we cannot exclude that dmrt1 and sox9 are expressed in different cell types. Stronger candidates for sox9 regulation in P. vitticeps are nr5a1 and fgf9, which are expressed highly in both sexes and could be responsible for the initial early and high sox9 expression in both sexes. The downregulation of sox9 in P. vitticeps female embryos at later stages is unlikely the result from decreases in nr5a1 and fgf9 expression in females, since their expression levels do not change at later stages. Another possibility is a female specific dominant repressor of sox9 expression. Two genes prominent in the female pathway in all vertebrates, foxl2 and cyp19a1, are expressed more strongly in females than males at stage 6/7, and the sex differential of both is increased at later stages: foxl2 by downregulation in male embryos at stage 12, and cyp19a1 by upregulation in females. The identity of a female specific dominant repressor of sox9 expression that acts to downregulate sox9 after the bipotential stage in female embryos remains speculative.
In mammalian systems, expression of AMH and its receptor AMHR initiates the regression of Müllerian ducts in male embryos . Post-natally AMH also has a role in controlling germ cell proliferation in both sexes [96, 97]. This latter role is seen also in fish and birds, and is thought to be the ancestral role of AMH in a vertebrate ancestor that lacked Müllerian ducts . However, in several fish, a copy of amh or amhr has a primary role in sex determination and defines a Y chromosome, while autosomal amh and amhr are involved in germ cell proliferation. Not surprisingly, the regulation of AMH is different in mammals and fish . In mammals it is induced by SOX9 activity, in combination with SF1 and GATA4 [99, 100]. In other systems, (chicken, T. scripta, and the American alligator Alligator mississippiensis) differential expression of amh was detected before differential expression of sox9 [2, 22, 101], although its action might be delayed by post-transcriptional regulation . P. vitticeps showed the non-mammalian pattern of early differential amh expression, which is high in both sexes but substantially higher in males. This suggests that in P. vitticeps, amh has a role in very early gonad differentiation as well as a later role in germ cell proliferation. P. vitticeps amh is not located on any of the four sex chromosome scaffolds identified in the P. vitticeps draft genome assembly . In fact, the scaffold with amh is not amongst those that have been cytogenetically mapped . An improved chromosome-scale genome assembly will determine whether or not amh is autosomal or encoded on the sex chromosomes.
The WNT/β-catenin signalling appears to be conserved among vertebrates as a key pathway in female gonad development that is repressed during male development, including fishes, reptiles with TSD, chickens and mammals. Commonly, WNT signalling activators RSPO1 and WNT4 finetune the strength of WNT signalling during sex differentiation being downregulated in males and/or upregulated in females [50, 103]. In P. vitticeps increased male expression of several WNT inhibitors is likely to result in the same outcome, lower protein levels of the effector β-catenin in males.
On the basis of our results, we propose a model for sex determination in P. vitticeps that is different to the mammalian system. In the classical model that was established in mammals, an upstream male determining factor initiates the runaway upregulation of SOX9 in males, triggering expression changes in key sex differentiation genes (including AMH). This does not apply to P. vitticeps, and probably not to reptiles generally. Early and central players for male development seem to be amh and dmrt1 rather than SRY and SOX9 as seen in mammals. Higher sox9 expression in P. vitticeps males is therefore a consequence, not a driver, of early sex differentiation. Nor does the model established for birds seem to apply, depending on dmrt1 and sox9 but not amh. Unlike in chicken, no known key male gene on the sex chromosomes showed a two-fold expression difference in gonads at the beginning of differentiation in P. vitticeps.
This and the high expression of male genes sox9 and amh in females after the onset of gonad differentiation led us to propose a dominance system for P. vitticeps in which the male trajectory is initially taken unless a dominant W-linked gene represses the male pathway to favour the female trajectory. In particular, expression of amh, dmrt1, WNT inhibitors, and foxl2 might be activated at a certain stage in gonad development (Fig. 5), regardless of the sex chromosome complement. GATA4 was shown to activate dmrt1  and a sex-independent activator was proposed to activate FOXL2 expression in mammals .
Which factor is responsible for the early and high expression of amh in P. vitticeps? In mammals it is SOX9, as well as SF1 and GATA4, amongst others [99, 100]. However, in some vertebrates amh expression preceded that of sox9 [2, 22, 101]. In our study both genes were highly expressed at the earliest examined stage, especially sox9, which subsequently displayed a dimorphic pattern. We cannot exclude the possibility that SOX9, together with GATA4 and SF1, acts as a sex non-specific activator of amh expression early in gonad development. However, we cannot exclude early SOX9 inhibition in the female pathway by preventing nuclear localisation via action of estrogen. In the absence of the W chromosome, in ZZ males, dmrt1 expression eventually supresses expression of foxl2 and the male pathway continues (Fig. 5); and outcompetes the female pathway. In ZW females, a dominant W-linked gene, directly or indirectly, would inhibit dmrt1 and WNT inhibitors, leading to establishment of WNT signalling and robust foxl2 expression which would tip the balance towards the female trajectory (Fig. 5).
Our study, the first RNA-sequencing of sex-typed GSD reptile embryos, provides a time series of embryonic gonad development in males and females. P. vitticeps shares downstream regulators of sex differentiation with reptile TSD species, and mammalian or avian GSD systems. However, those common regulators shift the timing of expression during gonad differentiation and the extent of sexually dimorphic expression. Expression profiles in P. vitticeps of key sex differentiation genes are more similar to known reptile TSD systems than they are to the mammalian GSD system. P. vitticeps has retained the central role for estrogen seen in lower vertebrates as high female specific expression of aromatase. However, P. vitticeps and possibly reptiles more broadly, seem to deviate from the mammalian expectation that early dimorphic upregulation of sox9 is integral to early male differentiation, highlighting the potential for plasticity and novelty in reptilian GSD. The candidate sex differentiation genes identified here, broaden our perspective on vertebrate sex differentiation by introducing a cast of new players and roles. This work is an essential step towards the identification of elusive master sex determining gene/s that govern sexual fate in reptiles with sex chromosomes.
Materials and methods
Animal breeding, egg incubations and embryo sampling
Animal breeding, egg incubation and embryo sampling was performed as described in  and was in accordance with animal ethics protocols of the University of Canberra and ARRIVE guidelines (https://arriveguidelines.org/). Here repeated, eggs were obtained during the 2017–18 breeding season from the research breeding colony at the University of Canberra. Breeding groups comprised three females (ZW) to two males (ZZ). Paternity was confirmed by SNP genotyping . Females were allowed to lay naturally, and eggs were collected at lay or within two hours of lay. Eggs were inspected for viability as indicated by presence of vasculature in the egg, and viable eggs were incubated at 28 °C in temperature-controlled incubators (± 1 °C) on damp vermiculite (4 parts water to 5 parts vermiculate by weight). Eggs were sampled at times corresponding to three developmental stages (6/7, 12 and 16). These stages equate to the beginning of gonad differentiation, recently differentiated gonad, and differentiated gonad, respectively, as described by Whiteley et al. [28, 29]. Between stages 4 to 8, gonads of both genotypes (ZW and ZZ) exhibit an elongated shape which develops into a rounder shape with defined cortex and medullary layers present at the end of this period and before gonad differentiation at stage 8/9. Testes differentiation is characterised by reduction of the cortex and proliferation of the medulla, within which seminiferous tubules form. Ovarian differentiation is characterised by a reduced medulla and a proliferating cortex with oogonia . Sample sizes are given in Fig. 1A. Embryos were euthanized by intracranial injection of 0.1 ml sodium pentobarbitone (60 mg/ml in isotonic saline). Individual gonads were dissected from the mesonephros under a dissection microscope and snap frozen in liquid nitrogen. Isolation of the gonad from the surrounding mesonephros was considered essential for studying transcriptional profiles within the gonad. Embryos were genotyped using previously established protocols [13, 28]. Briefly, this involved obtaining a blood sample from the vasculature on the inside of the eggshell on an FTA® Elute micro card (Whatman). DNA was extracted from the card following the manufacturer protocols, and PCR was used to amplify a W specific region  so allowing the identification of ZW and ZZ samples.
RNA extraction and sequencing
RNA extraction and sequencing was performed as described in . Here repeated, RNA from isolated gonad samples was extracted in randomized batches using the Qiagen RNeasy Micro Kit (Cat. No. 74004) according to the manufacturer protocols. RNA was eluted in 14 μl of RNase free water and frozen at -80 °C prior to sequencing. Sequencing libraries were prepared in randomized batches using 50 ng RNA input and the Roche NimbleGen KAPA Stranded mRNA-Seq Kit (Cat. No. KK8420). All sample RNA and library DNA was quantified using a Qubit Instrument (ThermoFisher Scientific, Scoresby, Australia), with fragment size and quality assessment using a Bioanalyzer (Agilent Technologies, Mulgrave, Australia). Nine randomly selected samples were sequenced per lane using the Illumina HiSeq 2500 system at the Kinghorn Centre for Clinical Genomics (Garvan Institute of Medical Research, Sydney) in a paired-end mode with a read length of 150 bp and 25 million read-pairs per sample were obtained on average.
For the gonadal transcriptome analysis of wild caught adult P. vitticeps ZZ and ZW individuals we used a previously published dataset . ENA accession numbers (https://www.ebi.ac.uk/ena) used for testes of ZZ male P. vitticeps: ERR753529, ERR413070; and for ovaries of ZW female P. vitticeps: ERR753530, ERR413082.
The paired-end RNA-seq libraries for P. vitticeps were trimmed using cutadapt  with -q 20 -m 20 -max-n 4 -trim-n. Trimmed reads were aligned to the P. vitticeps genome assembly pvi1.1 (GCA_900067755.1; http://ftp.ensembl.org/pub/release-100/fasta/pogona_vitticeps/dna/Pogona_vitticeps.pvi1.1.dna.toplevel.fa.gz)  using STAR (v2.7.0f)  with splice-aware alignment guided by the accompanying Ensembl gene annotation (pvi1.1.100). Parameters were chosen to output unique alignments (–outFilterMultimapNmax 1). To reduce the number of scaffolds in the genome assembly, scaffolds with no genes annotated were excluded.
Read count and differential gene expression analysis
Read count per gene was performed with htseq-count from the python package HTSeq  on the feature ‘exon’ (-t exon) with mode ‘union’ using the Ensembl gene annotation (pvi1.1.100) containing only genes tagged as protein-coding. Refer to Table S1 for raw counts per gene. In search for differentially expressed genes on the sex scaffolds, genes tagged as lncRNA (long non-coding RNA) were included as well (for Figure S3A). The R package DESeq2 (R 3.6.1, DESeq2 1.26.0 ) was used to perform differential gene expression analysis with default parameters. Genes with zero counts in all samples were excluded. Refer to Table S2 for complete results of a differential gene expression analysis between ZW and ZZ for each stage. To address gene expression changes from stage 6/7 to stage 12 differential gene expression analysis was performed between stages for each sex, ZZ and ZW. An adjusted p-value < 0.05 qualified the definition of differentially expressed gene (DEG) for a particular comparison. For normalisation the default method of DESeq2 was used, median of ratios, and whenever normalised read counts are plotted, this DESeq2 output was used. To compare expression levels of different genes the normalised read count was further normalised to the transcript length of the corresponding gene and expressed as nRPK (normalised read count per 1 kb transcript length). When more than one transcript was annotated for one gene, the average length of all annotated transcripts was used. If a gene of interest had no gene name assigned by the Ensembl gene annotation (pvi1.1.100), we looked for a gene name in the NCBI gene annotation (GCF_900067755.1 annotation release 100) by cross-referencing the Ensembl and NCBI gene IDs. Throughout the manuscript the Ensembl gene ID is given in addition to the gene name if a gene name was only assigned by the NCBI annotation and not by the Ensembl annotation. Such is the case for amh (ENSPVIG00000005953), col9a1 (ENSPVIG00000019420) and spink2 (ENSPVIG00000008173).
Histograms of gene expression height
The average nRPK of all samples per stage was used to plot histograms of gene expression. The histogram of stage 6/7 and stage 16 is very similar (Figure S10). Therefore, expression height of individual genes for stage 6/7 and stage 16 was plotted on the background of the stage 6/7 histogram. While the histogram for adults is more different (for example most genes in adults have an nRPK of roughly 300, whereas in embryonic stages it is roughly 100), individual gene expression height for adults was plotted on the background of the adult histogram.
Gene ontology analysis
Gene ontology enrichment analysis was done with ShinyGo v.0.66, http://bioinformatics.sdstate.edu/go/ , including inhere presented hierarchical clustering trees of gene ontology terms. As background served all protein-coding genes in the genome and the setting ‘Best matching species’ was used.
Weighted correlation network analysis
The R package WGCNA (R 4.0.3, WGCNA 1.70–3 ) was used. Read counts per gene (genes tagged protein-coding and lncRNA) were used as input (see above). Only embryonic stages were included. Only genes with a count > = 100 in at least two samples were included. For the gene network construction and identification of modules the following parameters were chosen: soft thresholding power = 10, type of network = unsigned (meaning that positive and negative correlations are considered), minimum module size = 20, merge cut height = 0.25. Figure panels S6D and E show the scale free topology model fit (signed R2) and the mean connectivity, respectively, over a range of soft thresholds.
For determination of the correlation between genotype and modules we have used the function ‘cor’ of the WGCNA package. To present positive correlations with the genotype ZZ, we have assigned ‘1’ for ZZ samples and ‘0’ for ZW samples (Figure S6B, left column). Vice versa, to present positive correlations with the genotype ZW, we have assigned ‘1’ for ZW samples and ‘0’ for ZZ samples (Figure S6B, right column). The two versions are mathematically redundant as they return same correlation values with opposite signs, and were chosen for presentational reason. Correlation values and p-values are given in Table S7.
For network visualisation Cytoscape (v. 3.9.1, ) was used. Visualizations are provided for modules that showed significant (p < 0.01) and high correlation (> 0.7) with ZZ or ZW genotype (Figures S7-9). The underlying nodes, edges and weights of those networks are provided in Table S7. Those modules were the 8/pink and 14/cyan having a strong correlation with ZW; and the 7/ black module having a strong correlation with ZZ. For meaningful visualisation only edges are represented with a weight > 0.08. Owing to the high number of nodes and edges in the black module, we provide a visualisation for edges with a weight > 0.11. The circular layout was chosen.
Gene nomenclature follows that of Kusumi et al. .
Availability of data and materials
Raw sequencing reads have been deposited in the Sequence Read Archive (SRA) repository of NCBI in BioProject PRJNA699086 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA699086).
Eggers S, Sinclair A. Mammalian sex determination—insights from humans and mice. Chromosom Res. 2012;20:215–38.
Hirst CE, Major AT, Smith CA. Sex determination and gonadal sex differentiation in the chicken model. Int J Dev Biol. 2018;62:153–66.
Ioannidis J, Taylor G, Zhao D, Liu L, Idoko-Akoh A, Gong D, Lovell-Badge R, Guioli S, McGrew MJ, Clinton M. Primary sex determination in birds depends on DMRT1 dosage, but gonadal sex does not determine adult secondary sex characteristics. Proc Natl Acad Sci. 2021;118: e2020909118.
Ayers KL, Davidson NM, Demiyah D, Roeszler KN, Grützner F, Sinclair AH, Oshlack A, Smith CA. RNA sequencing reveals sexually dimorphic gene expression before gonadal differentiation in chicken and allows comprehensive annotation of the W-chromosome. Genome Biol. 2013;14:R26.
Vining B, Ming Z, Bagheri-Fam S, Harley V. Diverse regulation but conserved function: SOX9 in vertebrate sex determination. Genes (Basel). 2021;12:486.
Nagahama Y, Chakraborty T, Paul-Prasanth B, Ohta K, Nakamura M. Sex determination, gonadal sex differentiation, and plasticity in vertebrate species. Physiol Rev. 2021;101:1237–308.
Sarre SD, Ezaz T, Georges A. Transitions between sex-determining systems in reptiles and amphibians. Annu Rev Genomics Hum Genet. 2011;12:391–406.
Valenzuela N. Evolution of the gene network underlying gonadogenesis in turtles with temperature-dependent and genotypic sex determination. Integr Comp Biol. 2008;48:476–85.
Greenbaum E, Carr JL. Sexual differentiation in the spiny softshell turtle (Apalone spinifera), a species with genetic sex determination. J Exp Zool. 2001;290:190–200.
Radhakrishnan S, Literman R, Neuwald J, Severin A, Valenzuela N. 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. PLoS ONE. 2017;12: e0172044.
Bista B, Wu Z, Literman R, Valenzuela N. Thermosensitive sex chromosome dosage compensation in ZZ/ZW softshell turtles, Apalone spinifera. Philos Trans R Soc B Biol Sci. 2021;376:20200101.
Literman R, Radhakrishnan S, Tamplin J, Burke R, Dresser C, Valenzuela N. Development of sexing primers in Glyptemys insculpta and Apalone spinifera turtles uncovers an XX/XY sex-determining system in the critically-endangered bog turtle Glyptemys muhlenbergii. Conserv Genet Resour. 2017;9:651–8.
Holleley CE, O’Meally D, Sarre SD, Marshall Graves JA, Ezaz T, Matsubara K, Azad B, Zhang X, Georges A. Sex reversal triggers the rapid transition from genetic to temperature-dependent sex. Nature. 2015;523:79–82.
Bull J, Vogt R. Temperature-dependent sex determination in turtles. Science. 1979;206:1186–8.
Capel B. Vertebrate sex determination: evolutionary plasticity of a fundamental switch. Nat Rev Genet. 2017;18:675–89.
Weber C, Zhou Y, Lee JG, Looger LL, Qian G, Ge C, Capel B. Temperature-dependent sex determination is mediated by pSTAT3 repression of Kdm6b. Science. 2020;368:303–6.
Ge C, Ye J, Weber C, Sun W, Zhang H, Zhou Y, Cai C, Qian G, Capel B. The histone demethylase KDM6B regulates temperature-dependent sex determination in a turtle species. Science. 2018;360:645–8.
Castelli MA, Whiteley SL, Georges A, Holleley CE. Cellular calcium and redox regulation: the mediator of vertebrate environmental sex determination? Biol Rev. 2020;95:680–95.
Whiteley SL, Holleley CE, Wagner S, Blackburn J, Deveson IW. JA Marshall Graves, A Georges, Two transcriptionally distinct pathways drive female development in a reptile with both genetic and temperature dependent sex determination. PLOS Genet. 2021;17:9465.
Yatsu R, Miyagawa S, Kohno S, Parrott BB, Yamaguchi K, Ogino Y, Miyakawa H, Lowers RH, Shigenobu S, Guillette LJ, Iguchi T. RNA-seq analysis of the gonadal transcriptome during Alligator mississippiensis temperature-dependent sex determination and differentiation. BMC Genomics. 2016;17:77.
Western PS, Harry JL, Marshall Graves JA, Sinclair AH. Temperature-dependent sex determination in the American alligator: expression of SF1, WT1 and DAX1 during gonadogenesis. Gene. 2000;241:223–32.
Czerwinski M, Natarajan A, Barske L, Looger LL, Capel B. A timecourse analysis of systemic and gonadal effects of temperature on sexual development of the red-eared slider turtle Trachemys scripta elegans. Dev Biol. 2016;420:166–77.
Lee L, Montiel EE, Navarro-Dominguez BM, Valenzuela N. Chromosomal rearrangements during turtle evolution altered the synteny of genes involved in vertebrate sex determination. Cytogenet Genome Res. 2019;157:77–88.
Deakin JE, Edwards MJ, Patel H, O’Meally D, Lian J, Stenhouse R, Ryan S, Livernois AM, Azad B, Holleley CE, Li Q, Georges A. Anchoring genome sequence to chromosomes of the central bearded dragon (Pogona vitticeps) enables reconstruction of ancestral squamate macrochromosomes and identifies sequence content of the Z chromosome. BMC Genomics. 2016;17:447.
Montiel EE, Badenhorst D, Tamplin J, Burke RL, Valenzuela N. Discovery of the youngest sex chromosomes reveals first case of convergent co-option of ancestral autosomes in turtles. Chromosoma. 2017;126:105–13.
Quinn AE, Georges A, Sarre SD, Guarino F, Ezaz T, Graves JAM. Temperature sex reversal implies sex gene dosage in a reptile. Science. 2007;316:411.
Ezaz T, Quinn AE, Miura I, Sarre SD, Georges A, Marshall Graves JA. The dragon lizard Pogona vitticeps has ZZ/ZW micro-sex chromosomes. Chromosom Res. 2005;13:763–76.
Whiteley SL, Holleley CE, Ruscoe WA, Castelli M, Whitehead DL, Lei J, Georges A, Weisbecker V. Sex determination mode does not affect body or genital development of the central bearded dragon (Pogona vitticeps). EvoDevo. 2017;8:25.
Whiteley SL, Weisbecker V, Georges A, Gauthier ARG, Whitehead DL, Holleley CE. Developmental asynchrony and antagonism of sex determination pathways in a lizard with temperature-induced sex reversal. Sci Rep. 2018;8:14892.
Frey N, Richardson JA, Olson EN. Calsarcins, a novel family of sarcomeric calcineurin-binding proteins. Proc Natl Acad Sci. 2000;97:14632–7.
Flynn R, Altier C. A macromolecular trafficking complex composed of β 2 -adrenergic receptors, A-Kinase Anchoring proteins and L-type calcium channels. J Recept Signal Transduct. 2013;33:172–6.
Marston S, Zamora JE. Troponin structure and function: a view of recent progress. J Muscle Res Cell Motil. 2020;41:71–89.
Zoli M, Pucci S, Vilella A, Gotti C. Neuronal and extraneuronal nicotinic acetylcholine receptors. Curr Neuropharmacol. 2018;16:338–49.
She Z-Y, Yang W-X. Molecular mechanisms involved in mammalian primary sex determination. J Mol Endocrinol. 2014;53:R21-37.
Hu Y-C, Okumura LM, Page DC. Gata4 is required for formation of the genital ridge in mice. PLoS Genet. 2013;9: e1003629.
Georges A, Li Q, Lian J, O’Meally D, Deakin J, Wang Z, Zhang P, Fujita M, Patel HR, Holleley CE, Zhou Y, Zhang X, Matsubara K, Waters P, Graves JAM, Sarre SD, Zhang G. High-coverage sequencing and annotated assembly of the genome of the Australian dragon lizard Pogona vitticeps. Gigascience. 2015;4:45.
Ionescu MI. Adenylate kinase: a ubiquitous enzyme correlated with medical conditions. Protein J. 2019;38:120–33.
Whiteley SL, Wagner S, Holleley CE, Deveson IW, Marshall Graves JA, Georges A. Truncated jarid2 and kdm6b transcripts are associated with temperature-induced sex reversal during development in a dragon lizard. Sci Adv. 2022;8:eabk0275.
Deveson IW, Holleley CE, Blackburn J, Marshall Graves JA, Mattick JS, Waters PD, Georges A. Differential intron retention in Jumonji chromatin modifier genes is implicated in reptile temperature-dependent sex determination. Sci Adv. 2017;3:e1700731.
Zhao Y, Mei Y, Chen HJ, Zhang LT, Wang H, Ji XS. Profiling expression changes of genes associated with temperature and sex during high temperature-induced masculinization in the Nile tilapia brain. Physiol Genomics. 2019;51:159–68.
Marroquín-Flores RA, Bowden RM, Paitz RT. Brief exposure to warm temperatures reduces intron retention in Kdm6b in a species with temperature-dependent sex determination. Biol Lett. 2021;17:20210167.
Kuroki S, Matoba S, Akiyoshi M, Matsumura Y, Miyachi H, Mise N, Abe K, Ogura A, Wilhelm D, Koopman P, Nozaki M, Kanai Y, Shinkai Y, Tachibana M. Epigenetic regulation of mouse sex determination by the histone demethylase Jmjd1a. Science. 2013;341:1106–9.
Omori Y, Kubo S, Kon T, Furuhashi M, Narita H, Kominami T, Ueno A, Tsutsumi R, Chaya T, Yamamoto H, Suetake I, Ueno S, Koseki H, Nakagawa A, Furukawa T. Samd7 is a cell type-specific PRC1 component essential for establishing retinal rod photoreceptor identity. Proc Natl Acad Sci. 2017;114:E8264–73.
Knower KC, Kelly S, Harley VR. Turning on the male – SRY, SOX9 and sex determination in mammals. Cytogenet Genome Res. 2003;101:185–98.
Colvin JS, Green RP, Schmahl J, Capel B, Ornitz DM. Male-to-female sex reversal in mice lacking fibroblast growth factor 9. Cell. 2001;104:875–89.
Bagheri-Fam S, Sim H, Bernard P, Jayakody I, Taketo MM, Scherer G, Harley VR. Loss of Fgfr2 leads to partial XY sex reversal. Dev Biol. 2008;314:71–83.
Rey R, Lukas-Croisier C, Lasala C, Bedecarrás P. AMH/MIS: what we know already about the gene, the protein and its regulation. Mol Cell Endocrinol. 2003;211:21–31.
Piprek RP. Genetic mechanisms underlying male sex determination in mammals. J Appl Genet. 2009;50:347–60.
Jeays-Ward K, Dandonneau M, Swain A. Wnt4 is required for proper male as well as female sexual development. Dev Biol. 2004;276:431–40.
Vainio S, Heikkilä M, Kispert A, Chin N, McMahon AP. Female development in mammals is regulated by Wnt-4 signalling. Nature. 1999;397:405–9.
Chassot A-A, Bradford ST, Auguste A, Gregoire EP, Pailhoux E, de Rooij DG, Schedl A, Chaboissier M-C. WNT4 and RSPO1 together are required for cell proliferation in the early mouse gonad. Development. 2012;139:4461–72.
Cruciat C-M, Niehrs C. Secreted and transmembrane Wnt inhibitors and activators. Cold Spring Harb Perspect Biol. 2013;5: a015081.
Bertho S, Pasquier J, Pan Q, Le Trionnaire G, Bobe J, Postlethwait JH, Pailhoux E, Schartl M, Herpin A, Guiguen Y. Foxl2 and its relatives are evolutionary conserved players in gonadal sex differentiation. Sex Dev. 2016;10:111–29.
Sánchez L, Chaouiya C. Primary sex determination of placental mammals: a modelling study uncovers dynamical developmental constraints in the formation of Sertoli and granulosa cells. BMC Syst Biol. 2016;10:37.
Cutting A, Chue J, Smith CA. Just how conserved is vertebrate sex determination? Dev Dyn. 2013;242:380–7.
Simpson ER, Clyne C, Rubin G, Boon WC, Robertson K, Britt K, Speed C, Jones M. Aromatase—a brief overview. Annu Rev Physiol. 2002;64:93–127.
Pask AJ. A role for estrogen in somatic cell fate of the mammalian gonad. Chromosom Res. 2012;20:239–45.
Antao AM, Ramakrishna S, Kim K-S. The role of Nkx3.1 in cancers and stemness. Int J Stem Cells. 2021;14:168–79.
A. Padmanabhan, V. Rao, A. M. De Marzo, C. J. Bieberich, Regulating NKX3.1 stability and function: post-translational modifications and structural determinants. Prostate. 2016; 76. https://doi.org/10.1002/pros.23144.
Arnesen C, Eich M-L, Pena MDCR, Cappel JR, Schwartz L, Rais-Bahrami S, Faraj SF, Prieto Granada C, Gordetsky JB. NKX3.1 and prostein expression in testicular tissue and sex cord-stromal tumors. Am J Surg Pathol. 2020;44:61–7.
Prescott JL, Blok L, Tindall DJ. Isolation and androgen regulation of the human homeobox cDNA, NKX3.1. Prostate. 1998;35:71–80.
Lee JM, Kim U, Yang H, Ryu B, Kim J, Sakuma T, Yamamoto T, Park J. TALEN-mediated generation of Nkx3.1 knockout rat model. Prostate. 2021;81:182–93.
Lee B, Park I, Jin S, Choi H, Kwon JT, Kim J, Jeong J, Cho B-N, Eddy EM, Cho C. Impaired spermatogenesis and fertility in mice carrying a mutation in the Spink2 gene expressed predominantly in testes. J Biol Chem. 2011;286:29108–17.
Piprek RP, Damulewicz M, Tassan J-P, Kloc M, Kubiak JZ. Transcriptome profiling reveals male- and female-specific gene expression pattern and novel gene candidates for the control of sex determination and gonad development in Xenopus laevis. Dev Genes Evol. 2019;229:53–72.
Perera EM, Martin H, Seeherunvong T, Kos L, Hughes IA, Hawkins JR, Berkovitz GD. Tescalcin, a novel gene encoding a putative EF-Hand Ca2+-Binding protein, Col9a3, and renin are expressed in the mouse testis during the early stages of gonadal differentiation**This work was supported in part by the Department of Pediatrics at the Univer. Endocrinology. 2001;142:455–63.
Piprek RP, Kolasa M, Podkowa D, Kloc M, Kubiak JZ. Transcriptional profiling validates involvement of extracellular matrix and proteinases genes in mouse gonad development. Mech Dev. 2018;149:9–19.
Sobinoff AP, Sutherland JM, Mclaughlin EA. Intracellular signalling during female gametogenesis. MHR Basic Sci Reprod Med. 2013;19:265–78.
Zup SL, Madden AMK. Gonadal hormone modulation of intracellular calcium as a mechanism of neuroprotection. Front Neuroendocrinol. 2016;42:40–52.
Aivatiadou E, Mattei E, Ceriani M, Tilia L, Berruti G. Impaired fertility and spermiogenetic disorders with loss of cell adhesion in male mice expressing an interfering Rap1 mutant. Mol Biol Cell. 2007;18:1530–42.
Bos JL, de Rooij J, Reedquist KA. Rap1 signalling: adhering to new models. Nat Rev Mol Cell Biol. 2001;2:369–77.
Matsumoto K, Ema M. Roles of VEGF-A signalling in development, regeneration, and tumours. J Biochem. 2014;156:1–10.
Singh N, Bhalla N. Moonlighting proteins. Annu Rev Genet. 2020;54:265–85.
Behrens J, von Kries JP, Kühl M, Bruhn L, Wedlich D, Grosschedl R, Birchmeier W. Functional interaction of β-catenin with the transcription factor LEF-1. Nature. 1996;382:638–42.
Zhao L, Wang C, Lehman ML, He M, An J, Svingen T, Spiller CM, Ng ET, Nelson CC, Koopman P. Transcriptomic analysis of mRNA expression and alternative splicing during mouse sex determination. Mol Cell Endocrinol. 2018;478:84–96.
Sreenivasan R, Jiang J, Wang X, Bártfai R, Kwan HY, Christoffels A, Orbán L. Gonad differentiation in Zebrafish is regulated by the canonical Wnt signaling pathway1. Biol Reprod. 2014;90:45.
Lawrence BM, O’Donnell L, Smith LB, Rebourcet D. New insights into testosterone biosynthesis: novel observations from HSD17B3 deficient mice. Int J Mol Sci. 2022;23:15555.
Heikkilä M, Prunskaite R, Naillat F, Itäranta P, Vuoristo J, Leppäluoto J, Peltoketo H, Vainio S. The partial female to male sex reversal in Wnt-4-Deficient females involves induced expression of testosterone biosynthetic genes and testosterone production, and depends on androgen action. Endocrinology. 2005;146:4016–23.
Roly ZY, Godini R, Estermann MA, Major AT, Pocock R, Smith CA. Transcriptional landscape of the embryonic chicken Müllerian duct. BMC Genomics. 2020;21:688.
Allbee AW, Rincon-Limas DE, Biteau B. Lmx1a is required for the development of the ovarian stem cell niche in Drosophila. Development. 2018;145:163394.
Tang B, Hu S, Ouyang Q, Wu T, Lu Y, Hu J, Hu B, Li L, Wang J. Comparative transcriptome analysis identifies crucial candidate genes and pathways in the hypothalamic-pituitary-gonadal axis during external genitalia development of male geese. BMC Genomics. 2022;23:136.
Kashimada K, Koopman P. Sry : the master switch in mammalian sex determination. Development. 2010;137:3921–30.
Hanley N, Hagan D, Clement-Jones M, Ball S, Strachan T, Salas-Cortés L, McElreavey K, Lindsay S, Robson S, Bullen P, Ostrer H, Wilson D. SRY, SOX9, and DAX1 expression patterns during human sex determination and gonadal development. Mech Dev. 2000;91:403–7.
Zhang X, Wagner S, Holleley CE, Deakin JE, Matsubara K, Deveson IW, O’Meally D, Patel HR, Ezaz T, Li Z, Wang C, Edwards M, Graves JAM, Georges A. Sex-specific splicing of Z- and W-borne nr5a1 alleles suggests sex determination is controlled by chromosome conformation. Proc Natl Acad Sci. 2022;119: e2116475119.
Munger SC, Aylor DL, Syed HA, Magwene PM, Threadgill DW, Capel B. Elucidation of the transcription network governing mammalian sex determination by exploiting strain-specific susceptibility to sex reversal. Genes Dev. 2009;23:2521–36.
Munger SC, Natarajan A, Looger LL, Ohler U, Capel B. Fine time course expression analysis identifies cascades of activation and repression and maps a putative regulator of mammalian sex determination. PLoS Genet. 2013;9: e1003630.
Shoemaker C, Ramsey M, Queen J, Crews D. Expression ofSox9, Mis, andDmrt1 in the gonad of a species with temperature-dependent sex determination. Dev Dyn. 2007;236:1055–63.
Barske LA, Capel B. Estrogen represses SOX9 during sex determination in the red-eared slider turtle Trachemys scripta. Dev Biol. 2010;341:305–14.
Nelson LR, Bulun SE. Estrogen production and action. J Am Acad Dermatol. 2001;45:S116-124.
Stewart MK, Mattiske DM, Pask AJ. Estrogen suppresses SOX9 and activates markers of female development in a human testis-derived cell line. BMC Mol Cell Biol. 2020;21:66.
Pask AJ, Calatayud NE, Shaw G, Wood WM, Renfree MB. Oestrogen blocks the nuclear entry of SOX9 in the developing gonad of a marsupial mammal. BMC Biol. 2010;8:113.
Malki S, Nef S, Notarnicola C, Thevenet L, Gasca S, Méjean C, Berta P, Poulat F, Boizet-Bonhoure B. Prostaglandin D2 induces nuclear import of the sex-determining factor SOX9 via its cAMP-PKA phosphorylation. EMBO J. 2005;24:1798–809.
Moniot B, Farhat A, Aritake K, Declosmenil F, Nef S, Eguchi N, Urade Y, Poulat F, Boizet-Bonhoure B. Hematopoietic prostaglandin D synthase (H-Pgds) is expressed in the early embryonic gonad and participates to the initial nuclear translocation of the SOX9 protein. Dev Dyn. 2011;240:2335–43.
Chen Y, Yu H, Pask AJ, Shaw G, Renfree MB. Prostaglandin D2 regulates SOX9 nuclear translocation during gonadal sex Determination in Tammar Wallaby. Macropus eugenii Sex Dev. 2017;11:143–50.
Zarkower D, Murphy MW. DMRT1: an ancient sexual regulator required for human gonadogenesis. Sex Dev. 2022;16:112–25.
Moses MM, Behringer RR. A gene regulatory network for Müllerian duct regression. Environ Epigenet. 2019;5:dvz017.
Durlinger ALL, Kramer P, Karels B, de Jong FH, Uilenbroek JTJ, Grootegoed JA, Themmen APN. Control of primordial follicle recruitment by anti-müllerian hormone in the mouse ovary. Endocrinology. 1999;140:5789–96.
Adolfi MC, Nakajima RT, Nóbrega RH, Schartl M. Intersex, hermaphroditism, and gonadal plasticity in vertebrates: evolution of the müllerian duct and Amh/Amhr2 signaling. Annu Rev Anim Biosci. 2019;7:149–72.
Mullen RD, Ontiveros AE, Moses MM, Behringer RR. AMH and AMHR2 mutations: a spectrum of reproductive phenotypes across vertebrate species. Dev Biol. 2019;455:1–9.
De Santa Barbara P, Bonneaud N, Boizet B, Desclozeaux M, Moniot B, Sudbeck P, Scherer G, Poulat F, Berta P. Direct Interaction of SRY-Related Protein SOX9 and steroidogenic factor 1 regulates transcription of the human anti-müllerian hormone gene. Mol Cell Biol. 1998;18:6653–65.
Tremblay JJ, Viger RS. Transcription factor GATA-4 enhances müllerian inhibiting substance gene transcription through a direct interaction with the nuclear receptor SF-1. Mol Endocrinol. 1999;13:1388–13401.
Western PS, Harry JL, Graves JAM, Sinclair AH. Temperature-dependent sex determination in the american alligator:AMH precedesSOX9 expression. Dev Dyn. 1999;216:411–9.
Lambeth LS, Ayers K, Cutting AD, Doran TJ, Sinclair AH, Smith CA. Anti-müllerian hormone is required for chicken embryonic urogenital system growth but not sexual differentiation. Biol Reprod. 2015;93:1–12.
Smith CA, Shoemaker CM, Roeszler KN, Queen J, Crews D, Sinclair AH. Cloning and expression of R-Spondin1in different vertebrates suggests a conserved role in ovarian development. BMC Dev Biol. 2008;8:72.
Lei N, Heckert LL. Gata4 regulates testis expression of Dmrt1. Mol Cell Biol. 2004;24:377–88.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17:10–2.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36:2628–9.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Kusumi K, Kulathinal RJ, Abzhanov A, Boissinot S, Crawford NG, Faircloth BC, Glenn TC, Janes DE, Losos JB, Menke DB, Poe S, Sanger TJ, Schneider CJ, Stapley J, Wade J, Wilson-Rawls J. Developing a community-based genetic nomenclature for anole lizards. BMC Genomics. 2011;12:554.
The animals used in this study were maintained and bred in the capable hands of Dr Wendy Ruscoe and Jacqui Richardson. We are grateful to the members of the wildlife genetics laboratory of the Institute for Applied Ecology for the many lively discussions on elements presented in this paper.
This work was supported by a Discovery Grant from the Australian Research Council (DP170101147) awarded to AG (lead), CEH, JD, JMG, Tariq Ezaz, Stephen Sarre, Lisa Schwanz and Paul Waters.
Ethics approval and consent to participate
All procedures were conducted with approval from the University of Canberra Animal Ethics Committee (AEC 17–17). We confirm that all the study was reported in accordance with ARRIVE guidelines (https://arriveguidelines.org/). All the methods were in accordance with animal ethics protocols of the University of Canberra.
Consent for publication
The authors declare no competing interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
FigureS1. Globaldifferential gene expression analysis in gonadogenesis. Figure S2.Muscle related genes are female specific in gonads at the beginning of differentiation.Figure S3. Sex determining gene not anymore dimorphic at stage 6/7 oryet has not been mapped to sex chromosomes. Figure S4. Early inhibitionof WNT signalling in male gonads by expression of WNT inhibitors. Figure S5.Dimorphic genes throughout all embryonic stages - candidates for novel sexdifferentiation genes. Figure S6. Weighted correlation networkexpression to identify female related gene sets. Figure S7. Networkvisualisation to reveal hub genes of female sex differentiation. Figure S8.Network visualisation to reveal hub genes of female sex differentiation. FigureS9. Network visualisation to reveal hub genes of male sex differentiation. FigureS10. Histograms of gene expression height per stage for comparison.
Raw counts per gene for all gonadal transcriptomes obtained with RNA-seq.
Results of a differential gene expression analysis of ZW versus ZZ for each stage.
Relevant genes in gonad development, especially transcription factors.
Expression of genes involved in WNT signalling.
Dimorphic expressed genes throughout all three examined embryonic stages.
Expression of genes involved in chromatin remodelling.
Modules of genes resulting from a weighted correlation network analysis.
About this article
Cite this article
Wagner, S., Whiteley, S.L., Castelli, M. et al. Gene expression of male pathway genes sox9 and amh during early sex differentiation in a reptile departs from the classical amniote model. BMC Genomics 24, 243 (2023). https://doi.org/10.1186/s12864-023-09334-0