Transcriptome sequencing and comparative analysis of adult ovary and testis identify potential gonadal maintenance-related genes in Mauremys reevesii with temperature-dependent sex determination

Mauremys reevesii is a classical organism with temperature-dependent sex determination (TSD). Gonad development in early life has recently received considerable attention but gonadal maintenance after sex differentiation in turtles with TSD remains a mystery. In this study, we sequenced the transcriptomes for the adult testis and ovary using RNA-seq, and 36,221 transcripts were identified. In total, 1,594 differentially expressed genes (DEGs) were identified where 756 DEGs were upregulated in the testis and 838 DEGs were upregulated in the ovary. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis suggested that the TGF-beta signaling pathway and Hedgehog signaling pathway have important roles in testis maintenance and spermatogenesis, whereas the Hippo signaling pathway and Wnt signaling pathway are likely to participate in ovary maintenance. We determined the existence of antagonistic networks containing significant specific-expressed genes and pathways related to gonadal maintenance and gametogenesis in the adult gonads of M. reevesii. The candidate gene Fibronectin type 3 and ankyrin repeat domains 1 (FANK1) might be involved with the regulation of testis spermatogenesis.


INTRODUCTION
Mauremys reevesii, also called Reeves' pond turtle, is one of the most common and widespread semiaquatic turtles in East Asia (Yin et al., 2016). This species is a classical organism with temperature-dependent sex determination (TSD). Previous studies developed a gene expression model for the embryonic gonads in reptiles with TSD (Yatsu et al., 2016;Czerwinski et al., 2016), and many of the sex-related genes have also been identified, for example, DMRT1 in Pelodiscus sinensis and Trachemys scripta Ge et al., 2017), and CYP19A1 in T. scripta (Matsumoto et al., 2013). In species with significant sex chromosomal differentiation, gonadal maintenance after sex differentiation depends on the expression of XX/XY or ZZ/ZW genes. However, species with TSD have almost identical genetic materials in the males and females. Thus, the mechanisms responsible for functional maintenance and gametogenesis during the sexual maturity period remain a mystery. Transcriptome sequencing provides a general representation of most of the transcripts that are expressed in specific cells or organs under particular conditions and stages (Liu et al., 2018), and it is considered the best option for identifying candidate genes in organisms that lack a reference genome (Yin et al., 2016).
In this study, we performed a comparative transcriptome analysis of the adult gonads (testis and ovary) in M. reevesii in order to identify molecular signaling cascades and gene expression networks. We also annotated all of the sequenced transcripts and enriched differentially expressed genes (DEGs), and pathways related to gonadal maintenance in these two sexual organs. These results provide crucial genomic information to facilitate further research into the regulatory mechanisms related to gonadal maintenance and gametogenesis during sexual maturity in turtle species with TSD.

Ethical approval
Procedures involving animals and their care were approved by the Animal Care and Use Committee of Anhui Normal University (#20170612).

Sample collection and RNA extraction
Six adult turtles (3♂, 3♀) were cultivated under the same rearing conditions (i.e., maintenance in a room at temperatures of 10-30 C) in the Provincial Key Laboratory for the Conservation and Exploitation Research of Biological Resources in Wuhu (31 33′N, 118 370′E, southeast of China) in 2017. The turtles were aged more than 6 years old in the sexual maturity stage, and the sexes of the samples were determined by dissection. Gonad samples (testicles and ovaries) were collected separately and dissected (Fig. 1). The gonad samples were flash frozen in liquid nitrogen and stored at -80 C until RNA was extracted (Yin et al., 2016).
RNA was extracted from each sample with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To effectively remove the genomic DNA, we added RNase-free DNase I (Takara, Dalian, China) to the reaction mixture for at least 10 min. The extracted RNA was quantified with the Nanodrop system (Thermo, Wilmington, DE, USA) and the fragment size distribution were checked by 1.5% agarose gel electrophoresis. RNA from each sample was used for cDNA synthesis and sequencing.

cDNA library construction and sequencing
We used 15 mg of total RNA for cDNA library construction using a TruSeq Ò RNA LT Sample Prep Kit v2 (Illumina, San Diego, CA, USA) according to the manufacturer's protocol.
Six cDNA libraries (3♂, 3♀) were constructed by performing end-repair, 3′-end adenylation, and adapter ligation and enrichment. Sequencing was performed using the Illumina Hiseq 2500 platform by Genergy Bio-technology Co. Ltd (Shanghai, China). All sequence reads have been deposited with the NCBI (GenBank accession No. SRP153785).
For samples with three biological replicates, the differential expression of unigenes between males and females was analyzed using the R (v3.4.2) package DESeq (Anders & Huber, 2010). Based on the negative binomial distribution model, DESeq provides statistical routines for determining differential expression in digital gene expression data (Hu et al., 2018). The statistical P-values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate (Benjamini & Hochberg, 1995). Unigenes with a corrected P-value < 0.05 and log 2 FoldChange ! 1 or log 2 FoldChange -1 were considered to be DEGs. Significant specific-expressed genes (SEGs) were identified from DEGs. Visualizations of the analyses, including the heatmap, volcano, GO, and KEGG enrichment results, were performed in the R package with ggplot2.

Validation of transcriptome data by real-time qRT-PCR
To verify the accuracy of the transcriptome data, three DEGs and nine SEGs were selected randomly for real-time qRT-PCR. All of the reactions were performed using three technical replicates and three biological replicates to validate the reliability of the results. A Real-time Detection System (C1000 Thermal Cycler; Bio-Rad, Hercules, CA, USA) was employed to perform real-time qRT-PCR. The PCR reaction system comprised the following: 10 mL SuperReal PreMix Plus (Tiangen, Beijing, China), 0.5 mL of each primer (10 mM in total), one mL of template cDNA, and eight mL of RNase-free ddH 2 O. The expression profile was detected in triplicate wells under the following protocol: 95 C for 15 min, 40 cycles for 10 s at 95 C and 20 s at 58 C, and 30 s at 72 C, before finally ramping from 65 to 95 C at 0.5 C per 5 s to generate a melting curve. The turtle GADPH gene was selected as the internal control gene. The relative expression level of each contig in different samples was calculated using the 2 -DDCT method (Livak & Schmittgen, 2001). The primers employed are listed in the supplemental materials (Table S1).

Sequencing and de novo assembly
Six libraries (3♂, 3♀) were established and sequenced (Table 1). In total, 508,785,280 paired-end raw reads with lengths over 100 bp were generated. We obtained 496,163,850 (97.57%) high-quality reads with an average length of 145.87 bp. The de novo assembly results yielded 557,690 isoforms, which clustered into 437,767 unigenes with an N50 of 631 bp. The GC content of the entire final assembly was 53.46%.

Unigene functional annotation
Among the total of 437,767 unigenes, 188,633 unigenes were annotated with the seven databases. The results are shown in Table S2. Among the databases, 20.08% of the unigenes were aligned to the Nr protein database with an e-value threshold of e-5 ( Fig. S1). In total, 17,361 GO annotations were assigned to three terms: biological processes (12,137 annotations), molecular functions (1,379 annotations), and cellular component (3,845 annotations). Among the biological processes, cellular, metabolic, organic substance metabolic, and primary metabolic processes were the most highly represented. The majority of the proteins assigned to molecular functions were associated with binding and catalytic activity. Among the cellular components, cell and intracellular proteins were the most highly represented (Fig. S2).
KEGG pathway analysis annotated 12,456 transcripts into 320 pathways, where 24.6% of all the annotated transcripts were related to metabolism, and 1,922 transcripts (15.4%) were involved with signaling molecules and interaction. The cytokine-cytokine receptor interaction pathway (ko04060; 804 sequences) was the most highly represented.

Analysis of DEGs
In total, 36,221 transcripts could be detected in the combined sex transcriptomes. After strict filtration, 756 DEGs were upregulated in the testis and 838 DEGs in the ovary (Fig. 2). The DEGs are listed in Table S3. A heatmap illustrating the hierarchical  the reliability and accuracy of the RNA-Seq method (Fig. 6). The results suggested that all of the genes had similar expression patterns according to both RNA-Seq and qRT-PCR, which also confirmed that RNA-Seq accurately quantified the expression of genes in the ovaries and testes.

DISCUSSION
In this study, we conducted the first transcriptomic gene expression analysis of adult gonads in a turtle species with TSD. In the gonads, we identify testis-and ovary-specific genes, and defined the downstream gene expression cascades for both the male and female pathways (Fig. 7). Existence of antagonistic signals and pathways related to gonadal maintenance and gametogenesis in M. reevesii In the sexual maturity period, steroid hormones such as estrogen and androgen are necessary for gonadal maintenance based on the sexually dimorphic expression of CYP19A1 and AR. CYP19A1, which is a gene encoding an enzyme that catalyzes conversion from androgens to estrogens (Piferrer, 2013;Matsumoto et al., 2013). CYP19A1 and other important genes such as SOX9 and FOXL2, which are related to sex determination and differences, exhibited sex-dimorphic expression in the early thermosensitive period (TSP) (Tang et al., 2017), and they still played core roles in the adult gonads in M. reevesii. We also identified conserved genes with known functions in sexual differentiation comprising AMH, GDF9, and WNT4. These genes also had obvious sexually dimorphic expression patterns. The male-specific expression of SOX9 mRNA during fetal and adult life indicates that it has a fundamental role in testis development in the turtle (Kent et al., 1996;She & Yang, 2017). In turtles, AMH is regulated by SOX9 and this pattern is consistent with SOX9 upregulating AMH in the mammalian testis (Bieser & Wibbels, 2014). Typical ovarian-specific typical markers (FOXL2 and CYP19A1) and testicular-specific marker SOX9 are often expressed in a mutually exclusive manner in the gonads (Matsumoto & Crews, 2012). FOXL2 suppresses SOX9 transcription by cooperatively binding with estrogen receptors in the regulatory region of SOX9 (Uhlenhaut et al., 2009). GDF9 can regulate a wide range of activities in the granulosa and theca cells, including the secretion of steroid hormones (Elvin, Yan & Matzuk, 2000), and it significantly suppresses the expression of AMH as well as being co-expressed with CYP19A1 during the sexual maturity period . WNT4 plays a key role in the female sexual development pathway by controlling steroidogenesis in the gonads and possibly supporting oocyte development (Vainio et al., 1999), as well as repressing typical male-specific steroidogenesis (Bernard & Harley, 2007). DMRT1 is a candidate master male sex-determining gene in T. scripta with TSD during early gonad development . However, in the adult turtle gonads, DMRT1 was detected at low levels and it did not exhibit transcriptional sexual dimorphism were significantly enriched in the testis (Fig. 7). The Hippo signaling pathway has critical roles in controlling cell proliferation, self-renewal, differentiation, and apoptosis in most tissues and organs in diverse species (Lyu et al., 2016) and it represents a molecular target for the regulation of mouse ovarian functional remodeling that could be used to regulate the proliferation and differentiation of ovarian function . Our analysis of M. reevesii suggested that Wnt signaling appears to function downstream of estrogen and the ovary-promoting effects of the Wnt signaling pathway may be functionally conserved in mammals and reptiles. Activation of the Wnt pathway antagonizes nuclear SOX9 expression (Mork & Capel, 2013). The TGF-beta signaling pathway plays an important role in the differentiation of male germ cells in non-mammalian vertebrates , for example, AMH is a growth factor in the TGF-beta family with a central role in testis formation. The Hedgehog signaling pathway has crucial roles in the development of diverse tissue and organ systems in the embryo, and in the regulation of adult tissues (Migone et al., 2017). Candidate gene FANK1 might be involved with the regulation of testis spermatogenesis FANK1 is exclusively expressed in the testis from the meiosis phase to the haploid phase of spermatogenesis in mice, and it may have a crucial functional role in spermatogenesis as a transcription factor (Dong et al., 2014). Indeed, reduced sperm numbers and increases in apoptotic spermatocytes were observed in FANK1 knockdown mice. In M. reevesii, FANK1 was strongly upregulated in the adult testis. Thus, it is possible that FANK1 play a pivotal role in spermatogenesis as a transcription factor in the adult gonads. Future studies should investigate the function of FANK1 during the regulation of testis maintenance.

CONCLUSIONS
According to our results, we propose a model where the genetic components related to steroid hormones comprise networks with multiple feedback loops. The networks are antagonistic in males and females, where differences in gene expression accumulate and converge to allow the antagonistic regulation of steroid hormones to maintain the appropriate balance in males and females. FANK1 could be involved with the regulation of testis spermatogenesis in the adult turtle gonads.