Introduction

The Chinese giant salamander Andrias davidianus is a typical urodele and an important species both as a biological resource and with respect to its value as a living fossil (Gao and Shubin 2003). The species was historically distributed in most areas of China, but environmental degradation has led to its severe decline in the wild, and it is classified as endangered by the International Union for Conservation of Nature and Nature Resources. Due to its palatability and nutritional value, A. davidianus is widely cultured for human consumption in many provinces of China. Since 2006, commercial sale of F2 generation Chinese giant salamander has been allowed by the Chinese government. In 2015, to expand the market, the requirement for a “Manage and Exploit License” and “Transport License” was discontinued and replaced with an identification card. Commercial breeding and rearing of giant salamander has become an important industry. However, sex growth dimorphism, with males growing more rapidly than females, has frustrated the industry to some extent. Sex differentiation is a fundamental biological process closely related to reproduction. To improve A. davidianus breeding technology, exploring the molecular mechanism underlying sex determination is critical. Amphibians played a key role as a transitional group in the evolution of vertebrates and possess one of largest and most diverse genome structures among vertebrates (Vinogradov 1998). Salamanders show a wide range of variation in genome size and karyotype. In Urodela, the largest amphibian genome and most complex karyotypes were exhibited in ancestral species (Voss et al. 2011). In salamander, the average genome is greater than 35 pg, with a high number of microchromosomes relative to other amphibians (Gregory 2002). Cryptobranchus alleganiensis possesses 30 chromosome pairs with both micro- and macro-chromosomes (Voss et al. 2011). The Ambystoma mexicanum genome consists of 14 chromosome pairs, and the estimated genome is ~ 32 Gb (Keinath et al. 2015). Andrias davidianus has 30 chromosome pairs including 19 pairs of microchromosomes (Zhu et al. 2002), with a genome of ~ 50 Gb (Geng et al. 2017). Transcriptome analysis is an efficient method of studying the large and complex genomes of amphibians. Using high-throughput sequencing, studies have been carried out on single and multiple genes. Transcriptome sequencing has been conducted in many fish species, including Cynoglossus semilaevis (Chen et al. 2014a, b), Pelteobagrus fulvidraco (Chen et al. 2014a, b), and Oreochromis niloticus (Tao et al. 2013). Due to the high cost of sequencing large genomes, transcriptome sequences are particularly useful for the study of amphibians. Transcriptome data of Anura and Urodela have been reported in species including Espadarana prosoblepon (Savage et al. 2014), Lithobates yavapaiensis (Savage et al. 2014), Lithobates clamitans (Robertson and Cornman 2014), Pseudacris regilla (Robertson and Cornman 2014), Rana chensinensis, Rana kukunoris (Yang et al. 2012), Xenopus (Silurana) tropicalis (Haselman et al. 2015), Xenopus laevis (Blowe et al. 2013), Hynobius chinensis (Che et al. 2014), Cynops pyrrhogaster (Nakamura et al. 2014), and Andrias davidianus (Li et al. 2015).

The sex determination mechanism of amphibians is complex. Sex determination involving male heterogamety (XX: XY) and female heterogamety (ZZ: ZW) has been found in both anurans and urodeles (Hillis and Green 1990). The W-linked DM domain gene (DMW), a Dmrt1 gene homolog, was determined to be the sex-determining gene (Yoshimoto et al. 2008) in X. laevis, but DMW has not been reported in X. (S.) tropicalis. In A. mexicanum, the major sex-determining factor and the single nucleotide polymorphism with alleles that associate with male and female phenotypes were localized to the genetic map (Smith et al. 2005, 2009). Sex-determining genes in other vertebrates include AMHY (Hattori et al. 2012), DmY (Matsuda et al. 2002), SRY (Koopman et al. 1991), FOXL2 (Boulanger et al. 2014), and Dmrt1 (Smith et al. 2009). The sex-determining gene is not highly conserved among vertebrates, and transcriptome sequences are useful in identifying candidate sex-determining genes.

Current studies of A. davidianus are primarily focused on immunity, with transcriptome information on several tissues available (Tao et al. 2013; Li et al. 2015; Huang et al. 2017). Several sex-differentiated genes have been reported (Sessions et al. 1982; Hu et al. 2016b).

The goal of the present study was to sequence the gonad transcriptome to identify ovary-biased and testis-biased genes and sex-related pathways and reveal the sex determination mechanism. Ten genes previously found to be involved in sex differentiation were isolated from the RNA-seq database, and RT-qPCR was performed to detect the expression profile in sex-reversed and normal individuals identified by a female-specific marker.

Materials and methods

Sample preparation

Ten (five males, five females) one-year-old apparently healthy A. davidianus specimens, were obtained from Zhejiang Yongqiang Chinese Giant Salamander Ltd. (Jinhua, Zhejiang Province, China). Phenotype sex was identified by histology. Ovaries and testes were dissected, and the samples were immersed in liquid nitrogen and stored at −80 °C for RNA extraction. This study was performed according to the guidelines of Yangtze River Fisheries Research Institute on the use of laboratory animals. Experimental protocols were approved by the Yangtze River Fisheries Research Institute Care Committee (No. 2013001).

RNA extraction, library construction, and sequencing

Total RNA was extracted using Trizol reagent (Invitrogen, USA) following manufacturer’s instructions and treated with RNase-free DNase I (Takara, China) for 30 min at 37 °C to remove residual DNA. RNA quality and quantity was tested using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and RNase-free agarose gel electrophoresis. Equal quantities of RNA from ovary and from testis of each salamander were, respectively, pooled for the transcriptome assembly library. Poly(A) mRNA was isolated using oligo-dT beads (Qiagen, Germany). Subsequently, the RNA was broken into short fragments, and the first-strand cDNA was generated by random hexamer-primed reverse transcription. Finally, the cDNA fragments were subjected to end reparation poly(A) addition and ligated to sequencing adapters. After purification and enrichment, the cDNA libraries were constructed and sequenced on the Illumina sequencing platform (Illumina HiSeq™ 2500).

De novo assembly and annotation

Before transcriptome assembly, all the raw reads were filtered to obtain clean reads. All reads containing adaptor sequences, reads with > 10% unknown bases, and sequences with > 50% of bases having a quality value ≤ 5 were removed. All the clean reads were assembled, and unigenes were detected using the Trinity Program (Grabherr et al. 2011). All unigenes were compared with the NCBI non-redundant (http://www.ncbi.nlm.nih.gov/), Swiss-Prot (http://www.expasy.ch/sprot), Clusters of Orthologous Groups of Proteins (COG) (http://www.ncbi.nlm.nih.gov/COG/), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, using the BLAST program with E value cut-off of 1e-5 (Altschul et al. 1997). Gene ontology (GO) annotation was analyzed by BLAST2GO software (Conesa et al. 2005), and the Web Gene Ontology (WEGO) tool was used to conduct function classification (Ye et al. 2006). The unigene sequences were scanned according to expression level. The threshold value was determined in multiple tests by the Benjamini-Yekutieli method (Benjamini and Yekutieli 2001). Criteria for differentially expressed genes (DEG) was a false discovery rate < 10−3 and expression level differing at least twofold between ovary and testis (Audic and Claverie 1997). Pathway enrichment investigation was carried out using the KEGG pathway database.

Characterization of the simple sequence repeat markers

Potential simple sequence repeat (SSR) markers were identified among unigenes by the MISA tool (http://pgrc.ipk-gatersleben.de/misa/). Simple sequence repeats were searched with motifs of sizes ranging from di- to hexa-nucleotides. Repeat units of SSRs were set as ≥ 6 repeat units for di-nucleotides, ≥ 5 repeat units for tri-nucelotides, and ≥ 4 repeat units for tetra-, penta-, and hexa-nucelotides. Separation distance of two SSRs > 100 bases were considered separate SSRs.

Identification of physiological sex and genetic sex

Genomic DNA was extracted from four females and double-digested using the restriction enzymes EcoRI and NlaIII (New England Biolabs USA). The RAD-seq libraries were constructed, and RAD-seq was performed as described by Recknagel et al. (2013). The female files were matched to the male genome, revealing 308 sequences that were not mapped to the male genome database. PCR amplification was carried out to validate the potential female-specific sequence. The potential marker Adf324 amplified a specific band only in females and was used to identify genetic females (Table S1). 80 individuals from normal group, 50 individuals from high temperature at 28 °C and 50 individuals from estradiol exposure were killed and the gonads were collected and fixed in 4% paraformaldehyde (pH 7.5) for 24 h, and subsequently stored in 70% ethanol. Gonads were dehydrated in an ethanol series, cleared in xylene, and embedded in paraffin. Sections were cut at 5–6 µm and stained with hematoxylin–eosin. The physiological sex was identified under light microscopy (Olympus). The genetic sex was identified by the female-specific marker (Adf324) previously identified (Table S1). The PCR reaction was conducted on 20 µL using 1.1XT3 Super PCR Mix (TsingKe, Beijing, China) according to the manufacturer’s instructions. The PCR reactions were performed under the following conditions: 94 °C for 3 min, 35 cycles of 94 °C for 30 s, 59.5 °C for 30 s, and 72 °C for 40 s; then 72 °C for 5 min. When physiological sex was male and genetic sex female, the specimen was considered neomale, and when physiological sex was female and genetic sex male, it was considered neofemale.

RT-qPCR analysis

From the differentially expressed genes, we selected 27 well-known candidate genes involved in sex differentiation or gonad development from the comparative transcriptome. We randomly selected 10 of the 27 genes and evaluated the expression profile in gonad of female, male, neomale, and neofemale by RT-qPCR. Primers were designed according to the transcriptome sequences (Table S1). cDNA was synthesized using a QuantScript RT kit (Takara, Dalian, China) with RNA extracted from gonads of an eighteen-month-old A. davidianus according to manufacturer’s protocol. Quantitative real-time PCR was performed using the ABI Q5 real-time PCR system (Applied Biosystems, California, USA). The amplification protocol was as previously described (Hu et al. 2016a), and β-actin was used as the internal control.

Results

Transcriptome sequencing and de novo assembly

cDNA libraries were constructed from mRNA of ovary and testis. Transcriptome sequencing was performed on the Illumina sequencing platform. Raw reads were deposited in the GenBank database under numbers SRR3308418 and SRR3308420. A total of 40,783,222 and 46,128,902 raw reads were obtained with 5,097,902,750 and 5,766,112,750 bp in ovary and testis, respectively (Table 1). All clean reads were assembled by the Trinity program without the reference genome sequence. Approximately, 70% of the reads could be mapped back to the assembled transcripts. The N50 length of the unigenes was 1557 bp, and the maximum and minimum lengths of the unigenes were 14,449 bp and 201 bp, respectively (Table S2). There were 60,825 (75.56%) unigenes ranging from 200 to 1000 bp; 11,075 (13.76%) from 1,000 to 2,000 bp; 4,774 (5.93%) from 2000 to 3000 bp; and 3823 (4.75%) unigenes of more than 3000 bp (Fig. 1).

Table 1 Number of original and clean reads mapped to their respective assembly transcript of Andrias davidianus
Full size table
Fig. 1
figure 1

Length distribution of assembled unigenes

Full size image

Functional annotation and classification of the Andrias davidianus transcriptome

Using the BLAST program (E value < 1e−5), the NR, Swissprot, COG, and KEGG databases were employed to annotate all assembled unigenes. A total of 36,211 and 30,789 unigenes were classified in the NR and Swissprot databases, respectively, while 11,075 and 18,318 unigenes were assigned to the COG and KEGG databases, respectively. A total of 36,609 (45.48% of 80,497) unigenes had best hit to one or more of the databases. The Gene Ontology (GO) project defines three categories: biological processes, molecular function, and cellular components. According to the sequence homology, 18,971 (23.56% of 80,497) unigenes were assigned to GO categories (Fig. 2). Cellular process (12,483; 65.81% of 18,971), metabolic process (10,414; 54.89%), and biological regulation (5782; 30.48%) accounted for the largest segment of the biological processes category. Cell (13,490; 71.11%), cell part (13,489, 71.10%), and organelle (9680; 51.03%) were the most common cellular components. Binding (11,568; 60.98%) and catalytic activity (7,465; 39.35%) dominated in molecular functions. All unigenes were aligned to the COG database. Results assigned 11,075 (13.75% of 80,497) unigenes to a COG classification, grouped into 24 function categories (Fig. 3). The largest category was general function prediction only (3,437; 18.31%), followed by post-translational modification, protein turnover, and chaperones (1,502; 8.00%); and replication, recombination, and repair (1,254; 6.61%).

Fig. 2
figure 2

Gene ontology (GO) classification of Andrias davidianus

Full size image
Fig. 3
figure 3

Cluster of orthologous groups (COG) functional classification of Andrias davidianus

Full size image

Sex-biased genes and sexual dimorphic pathways

A total of 15,615 female-biased and 19,679 male-biased unigenes were indentified from the gonad transcriptome, showing significantly different expression in male and female (Table S3). From known genes, 8907 (57.04%) female-biased and 10,385 (52.77%) male-biased unigenes were annotated. In addition, 3669 unigenes were exclusively expressed in male and 1830 unigenes were exclusively expressed in female (Table S3). From the sex-biased genes, we searched 27 well-known candidate genes characterized as sexually dimorphic in other species, of which nine were female-biased and 18 were male-biased (Table 2). To further elucidate the biological pathways involved in sex-biased genes, 10,101 (28.62%) of 35,294 DEGs assigned to the 240 KEGG pathways were analyzed. We found a greater number of sex-biased genes, including 1,571 unigenes, in the metabolic pathway, compared to the other pathways.

Table 2 Differentially expressed genes related to sex determination and differentiation in Andrias davidianus
Full size table

Several pathways were involved in gonad development, including the Wnt signaling pathway (Fig. 4), progesterone-mediated oocyte maturation, GnRH signaling pathway, TGF-beta signaling pathway, steroid hormone biosynthesis, and steroid biosynthesis. In the Wnt signaling pathways, 63 genes were identified, including 54 DEGs comprising 20 female-biased, 12 male-biased, and 22 sex-biased genes (Fig. 4).

Fig. 4
figure 4

Wnt signaling pathway. The purple background indicates genes or pathways identified in A. davidianus; red background indicates ovary-biased genes or pathways; green background indicates testis-biased genes or pathways; blue background indicates pathways containing ovary-biased and testis-biased genes. (Color figure online)

Full size image

Characterization of simple sequence repeat (SSR) markers

To develop an SSR marker for A. davidianus using MISA software, we identified 3,427 SSRs from 80,497 examined sequences, with 3,111 containing SSR sequences, which accounted for 3.86% of the total examined sequences (Table 3). Multiple SSRs were observed in 271 sequences, and 146 SSRs were present in compound formation. Five types of SSRs were detected, with the di-nucleotide repeat motif comprising the largest group (1,852, 54.04%), followed by the tri-nucleotide (1,155, 33.70%), tetra-nucleotide (350, 10.21%), penta-nucleotide (40, 1.17%), and hexa-nucleotide (30, 0.88%) (Fig. S1).

Table 3 SSR summary from the Andrias davidianus transcriptome
Full size table

Identification of physiological sex and genetic sex

Physiological sex was identified by histology sections. When follicles were present, the individual was considered female (Fig. 5a). When testicular lobule and spermatogonia were observed, the individual was identified as male (Fig. 5b). 24 female and 24 male individuals from normal group, 12 female and 12 male individuals from high temperature group, and 20 female and four male individuals from estradiol-exposed group were selected for the genetic sex identification. The results of the genetic sex showed that twenty-four females and 24 males were distinguished by presence or absence of a female-specific marker (Fig. 5c). Female-specific bands were observed in the 12 females and in three of 12 males exposed to 28 °C (Fig. 5d). The physiological sex of the males showing the female-specific band was different from the genetic sex, indicating sex reversal from female to male. In the estradiol-exposed groups, female-specific bands were detected in seven specimens among the 20 females, with none in the four males (Fig. 5e). Thirteen female specimens did not exhibit the female-specific band, suggesting reversal from male to female.

Fig. 5
figure 5

Identification of physiological sex and genetic sex. a Histology sections showing the female gonad. b Histology sections showing the male gonad. c Genetic sex was identified by female-specific marker, and 24 physiological females showed the clear female-specific band, while 24 physiological males did not. d Genetic sex of individuals exposed to high temperature (28 °C) were identified. e Genetic sex of individuals exposed to estradiol were identified. GC germ cell, GrC granulosa cells, SL seminiferous lobule, Sp spermatogonia, FC follicular cavity, SeC sertoli cell

Full size image

Expression of differentiated expressed gene

RT-qPCR revealed the expression profile of ten genes in gonad of male and female in accordance with the results of transcriptome sequencing. Male-biased genes Amh, Dmrt1, Fgf17, Gata4, and Sf1 exhibited up-regulation in gonad of female to male reversal, with no significant expression observed in gonad of female and male to female reversal (Fig. 6). Female-biased genes Zp3, CYP19a, CYP26a, TfIII, and Wnt2b showed down-regulation in gonad of male and reversed male. Expression of female-biased genes in gonad of female and reversed female showed no significant differences, except for TfIIIa, which was significantly down-regulated in reversed female (Fig. 6).

Fig. 6
figure 6

Relative expression profiles of DEG in gonad. F female, M male, RF reversed female, RM reversed male

Full size image

Discussion

Because of their transitional evolutionary status, the sex-determining mechanism of amphibians, including genetic sex determination (GSD) and environmental sex determination (ESD) has attracted the interest of researchers (Wallace et al. 1999). In A. davidianus, the sex ratio was 1:1 at 20 °C. Masculinization was induced in A. davidianus exposed to high temperatures during the sex differentiation stage (Hu et al. unpublished data), suggesting that A. davidianus possesses both GSD and ESD mechanisms. In previous studies, several sex-related genes were cloned, and their expression was analyzed (Hu et al. 2015; Hu et al. 2016b). However, due to species differences in sex-determining gene, transcriptome sequence screening of DEGs in males and females is more informative.

The N50 length was used to evaluate the assembly quality. The N50 length of our assembly was higher than observed in most other amphibians, including R. chensinensis (1333 bp) (Yang et al. 2012), H. arborea (700 bp) (Brelsford et al. 2013), R. kukunoris (1485 bp) (Yang et al. 2012), L. clamitan (686 bp) (Robertson and Cornman 2014) and P. regilla (687 bp) (Robertson and Cornman 2014). The annotation rate of A. davidianus gonad transcriptomes was higher than reported in C. alboguttata (22,695 of 68,947) (Reilly et al. 2011)and lower than in O. margaretae (18,933 of 37,906) (Tao et al. 2013).

In the present study, 15,615 female-biased and 19,679 male-biased unigenes in gonad transcriptomes were scanned (Table S3). Nearly 57.04% of the female-biased and 52.77% of the male-biased unigenes were annotated. In annotated male-biased genes, well-documented male genes were identified, including Dmrt1, AMH, Gata4, SF1, and Sox5. In chicken, reduction of dmrt1 protein expression by RNA interference was reported to lead to feminization of the gonads in male embryos. Affected males showed partial sex reversal and a decline of the testicular marker gene, while the ovarian marker gene was activated (Smith et al. 2009). In Odontesthes hatcheri, a duplicated copy of the anti-Müllerian hormone (amhy) is expressed much earlier in gonad development than is the amh gene. Following amhy gene knockdown in XY males, foxl2 and CYP19a1a mRNA expression was up-regulated, accompanied by development of ovaries (Hattori et al. 2012). In X. laevis, the DMW gene, paralogue of Dmrt1, is expressed exclusively in primordial gonads of ZW tadpoles (Yoshimoto et al. 2008). In the early sex determination stage, DMW mRNA expression is reported to be higher than Dmrt1. Transgenic tadpoles with the DMW expression vector showed ovarian cavities and primary oocytes (Yoshimoto et al. 2008). In DMW knockdown X. laevis, the gonad of the ZW female has been reported to show regressed germ cells, underdeveloped ovaries, or no oocytes (Yoshimoto et al. 2008). Hence, the DMW gene is considered the ovary-determining gene in X. laevis. Fushi tarazu factor-1 (FTZ-F1) is a nuclear receptor and an essential regulator of sex determination in mammals via action on the expression of steroidogenic enzymes (Val et al. 2003). In Oreochromis niloticus, Ftz-f1 is present in the 5′ flanking region of CYP19a1, suggesting Ftz-f1 as possibly engaged in regulation of CYP19a1 expression (Yoshiura et al. 2003). Wang et al. (2007) reported FOXL2 to interact with the ligand binding domain of Ad4BP/SF-1 to form a heterodimer and enhance Ad4BP/SF-1-mediated CYP19a1 transcription. However, SF1 has shown no effect on the expression of CYP19 in HEK293 and A6 cells (Oshima et al. 2006). In previous studies, we detected higher CYP19a expression in A. davidianus ovary than in testis, while SF1 exhibited a different expression pattern from CYP19a in gonad (Hu et al. 2016a, b), suggesting that SF1 does not up-regulate expression of CYP19a in A. davidianus. A similar expression profile was found in Trachemys scripta. CYP19 expression increases during ovarian development, while SF1 is more highly expressed during testis development (Ramsey et al. 2007). Feminizing temperatures and exogenous estrogen have been shown to inhibit SF1 expression, whereas inhibition of CYP19 and male-producing temperatures increase SF1 expression (Ramsey et al. 2007).

We scanned several documented female genes including CYP19a1, Zp3, Zp4, CYP26a. CYP19a1 produces aromatase, which irreversibly converts androgens to estrogens. In humans, the insertion of a mutation into the fifth exon of CYP19a1 resulted in hypoplastic ovaries, suggesting that CYP19a1 is necessary for ovary development (Akçurin et al. 2016). In R. rugosa, CYP19a expression was up-regulated in the undifferentiated gonad of female tadpoles during sex determination (Maruo et al. 2008). Expression of CYP19 was reported to be associated with sexual differentiation in Danio rerio (Trant et al. 2001). A recent study showed the CYP19a1a gene to perform functions in both male and female sex differentiation (Yin et al. 2017). In vertebrates, sex is differentiated by the timing of meiosis initiation. In mammals, germ cells in ovary enter meiosis earlier than germ cells in testis, while in Cyp26b1-knockout mouse embryo testis, germ cells enter meiosis precociously (Bowles et al. 2006). In Oreochromis niloticus, treatment with retinoic acid synthase inhibitor delayed meiotic initiation and resulted in down-regulation of CYP26a1 (Feng et al. 2015). CYP26a expression was induced in Acanthopagrus schlegelii treated with estradiol (Lau et al. 2013). The zona pellucida is essential for oogenesis and fertilization in mammals. Mutations of Zp2 and Zp3 can cause female infertility in humans (Liu et al. 2017). In Mus musculus, knockdown of Zp3 dramatically inhibited germinal vesicle breakdown and reduced the percentage of MII oocytes (Gao et al. 2017). Gonads are mainly responsible for producing gametes and sex hormones. Sex hormones bind to nuclear receptors to take part in pathways for specific physiological responses, including sex differentiation and germ cell development.

In the present study, we detected numerous sex-biased genes involved in the WNT signaling pathway. Two of 19 members of the Wnt family have been identified in follicular granulosa cells. Expression of Wnt2 has been detected in all stages of follicular development, but Wnt2-null mice exhibited no fertility defects, suggesting functional redundancy of Wnt2 (Monkley et al. 1996). Expression of Wnt4 was detected in embryogenesis and follicle formation (Vainio et al. 1999a, b), and studies of Wnt4 null mice showed it to be essential for female fertility (Boerboom 2009). Loss or gain of function of important components of the WNT pathway such as Wnt4, Fzd4, and Ctnnb1 results in abnormal folliculogenesis and steroidogenesis (Boyer et al. 2010). The Wnt signaling pathway was highly conserved during evolution and plays a vital role in gonad development. In mammals, the WNT signaling pathway has been shown to be critical to female development (Vainio et al. 1999a, b). In the red-eared slider turtle Trachemys scripta elegans, ectopic activation of the WNT signaling pathway in male gonads induces sex reversal from male to female, while inhibition of the pathway in ovary produced reversal from female to male (Mork and Capel 2013). Studies on fish have shown that Wnt/β-catenin signaling plays a crucial role in ovary differentiation in rainbow trout Oncorhynchus mykiss and the protandrous black porgy Acanthopagrus schlegelii (Wu and Chang 2009; Nicol and Guiguen 2011).

Temperature and sex hormone influence on amphibian sex differentiation has been demonstrated in many species. In Pleurodeles waltl, ZW neomales were produced by rearing larvae at elevated temperatures (Dournon and Houillon 1985). Estrogen was shown to induce feminization in urodele and anuran larvae, while androgen had no effect, or produced females, in urodeles (Alam et al. 2008). In this study, in 12 high temperature-treated males, we detected three genetic females (neomales). In 20 estrogen-treated females, we detected seven genetic females and thirteen neofemales. Expression of selected DEGs was detected in gonads of males, females, neomales, and neofemales. Male-biased genes were up-regulated in males and neomales, while no significant expression was observed in females and neofemales. The opposite expression profile was observed in female-biased genes (Fig. 6).

Similar results have been observed in fish. In Epinephelus merra, dmrt1 expression was up-regulated, and FOXL2 expression was down-regulated, during female-to-male sex reversal, and no significant expression difference was observed between male and neomale or female and neofemale (Alam et al. 2008). In Cynoglossus semilaevis, SOX9 expression was high in males and neomales exposed to elevated temperature and testosterone, while low expression was observed in females and neofemales. The female-biased gene FOXL2 exhibited the opposite expression profile (Dong et al. 2011). In chicken, expression of dmrt1 was up-regulated during female to male sex reversal (Smith et al. 2003). Results suggested that DEGs play a key role in determining and maintaining sex differentiation. Further study is required for confirmation.

We screened genes critical to sex determination in Andrias davidianus, for the first time, by comparative transcriptome sequencing of ovary and testis. We detected numerous sex-biased genes involved in pathways related to gonad development, including WNT signaling pathways. The female-specific marker was used to identify sex reversed and normal individuals. Expression profiles of ten sex-biased genes were validated in sex-reversed and normal individuals. These results expand understanding of the regulatory network associated with sex differentiation in A. davidianus, which will produce benefits for its artificial culture.