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

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 μg 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).

Sequence data processing and de novo assembly

The raw reads generated by Illumina sequencing were subjected to a screening process where the adapter sequences and low-quality reads were trimmed from the raw reads using Trim Galore v0.5.0 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) software. FastQC v0.11.8 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to check the quality of the pretreated data. Trinity v2.8.4 (https://github.com/trinityrnaseq/trinityrnaseq/releases/tag/Trinity-v2.8.4) was used to perform de novo assembly with the default parameters (Grabherr et al., 2011). The clean high quality data were spliced to obtain the reference sequences (transcripts) for subsequent analyses.

Gene annotation

Multiple public databases were used for homology annotation, including non-redundant nucleotide and protein sequences (NT/NR) (http://www.ncbi.nlm.nih.gov), Swiss-Prot (http://www.expasy.ch/sprot), euKaryotic Orthologous Groups (KOG) (https://ftp.ncbi.nih.gov/pub/COG/KOG/) based on BLASTX (NCBI blast 2.2.28) with an e-value of 1e-5 (Altschul et al., 1990), Gene Ontology (GO) (http://www.geneontology.org/) based on Blast2GO v2.5 with an e-value of 1e-6, and Kyoto encyclopedia of genes and genomes (KEGG; http://www.genome.jp/kegg) based on KAAS with an e-value of 1e-10 (Kanehisa et al., 2008). InterProScan v5.31-70.0 (http://www.ebi.ac.uk/interpro/download.html) was also used to predict the functional domains, signal peptides, and other protein characteristics by BLASTing against the Conserved Domain Database Interpro under the default parameters (Zdobnov & Apweiler, 2001). GO terms were assigned to three categories (biological process, molecular function, and cellular component) (Conesa et al., 2005). The KEGG Automatic Annotation Server was employed for KEGG annotation (Moriya et al., 2007).

DEG analysis

Clean reads were aligned using Bowtie 2 v2.3.4.3 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) (Langmead, 2010). Transcript abundance was estimated using RSEM v1.3.1 (http://deweylab.github.io/RSEM/) in the Trinity package (Li & Dewey, 2011). The expression profiles were detected for the transcripts using the gene expression method as the fragments per kilobase of transcript per million mapped reads (Trapnell et al., 2010). 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₂FoldChange ≥ 1 or log₂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 μL SuperReal PreMix Plus (Tiangen, Beijing, China), 0.5 μL of each primer (10 μM in total), one μL of template cDNA, and eight μL of RNase-free ddH₂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^−ΔΔCT 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 clustering of the DEGs was generated to visualize the overall gene expression pattern (Fig. 3). Among the 1,594 DEGs, we found nine significant SEGs, as shown in Table 2.

Enrichment analysis of DEGs

The 1,175 GO terms were annotated as biological processes (859; 73.1%), cellular component (109; 9.3%), and molecular function (207; 17.6%). DEGs were assigned to 39 processes by KEGG enrichment analysis, including Cellular Processes (310 DEGs, 22.27%), Metabolism (88 DEGs, 22.27%), Genetic Information Processing (216 DEGs, 22.27%), Environmental Information Processing (352 DEGs, 22.27%), Human Diseases (662 DEGs, 22.27%), and Organismal Systems (288 DEGs, 22.27%). We used the same strategies (Du et al., 2017) used to identify genes and pathways related to gonadal development and the gametogenesis system, and 20 significantly enriched GO terms were detected in the ovary and testis. These enriched sex-related terms included male gonad development, tube development in testis, embryo development, iron ion binding, and regulation of signal transduction in ovary. The 20 GO terms are shown in Fig. 4. In addition, 10 KEGG terms were significantly enriched in the ovary and testis, that is, the Hedgehog signaling pathway, pathways in cancer, TGF-beta signaling pathway, prostate cancer in testis, metabolic pathways, Hippo signaling pathway, ovarian steroidogenesis, steroid hormone biosynthesis, Wnt signaling pathway, and oocyte meiosis in ovary. The 10 KEGG terms are shown in Fig. 5.

Validation of transcriptome data by real-time qRT-PCR

A total of 12 genes comprising six ovary SEGs (FGF9, FOXL2, GDF9, WNT4, CYP19A1, and SAMD9), three testis SEGs (SOX9, AMH, and AR), and three DEGS (ARHGEF9, CKAP5, and PTH2R) were selected for validation using qRT-PCR in order to confirm 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.

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 (Wang et al., 2017). 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 (Ge et al., 2017). However, in the adult turtle gonads, DMRT1 was detected at low levels and it did not exhibit transcriptional sexual dimorphism between the testis and ovary. Our results suggest that DMRT1 does no play an important role in the adult differentiated gonads. Czerwinski et al. (2016) identified six male-producing temperature-enriched genes and four female-producing temperature-enriched genes with differences in expression prior to DMRT1 and CYP19A1 in the early TSP. The detection of these 10 genes suggests that a gonad-specific regulatory mechanism exists but it is not clear if it is causal. In our study, NOV, KDM6B, RBM20, PCSK6, AVIL, and TWIST1 were detected at low levels and their expression levels did not differ between the adult testis and ovary. HSPB6 was not detected. VWA2 and FDXR were upregulated in the ovary, but their expression levels were low in the testis, and only FANK1 was strongly upregulated in the adult testis. Our results suggest that the potential gonad-specific regulatory mechanism might differ between the embryo and adult.

The Hippo signaling pathway and Wnt signaling pathway were significantly enriched in the ovary, whereas the TGF-beta signaling pathway and Hedgehog signaling pathway 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 (Ye et al., 2017). 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 (Zhang et al., 2016), 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.