Ethics statement, experimental fish, and sample collection.

All procedures involving the handling and treatment of fish used during this study were approved by the Animal Care and Use Committee at Heilongjiang River Fisheries Research Institute. The fish used for this project were 20 male yellow catfish (10 1-yr-old juveniles, 10 2-yr-old adults) and 20 female yellow catfish (10 1-yr-old juveniles, 10 2-yr-old adults). All 40 fish came from one full-sibling family located at National Fish Original Species Farm in Zhaodong city, Heilongjiang Province, China. One-year-old fish and two-year-old fish were collected on August 16, 2012 and August 14, 2013, respectively. All fish sex was confirmed anatomically. The experimental fish were euthanized with tricaine methanesulfonate at 250 mg/l before sample collection. Gonad tissues from two different ages were collected and stored in 1.5 ml RNAlater tubes (Qiagen, Hilden, Germany) and then transferred to a −80°C freezer until prior to RNA extraction.

RNA extraction.

Samples were removed from the −80°C freezer and homogenized with TissueRuptor (Qiagen) to a fine solution. Total RNA were extracted from each sample with the RNeasy Lipid Tissue Kit (Qiagen) according to the manufacturer's instructions after the RNase free DNase I (Qiagen) treatment of genomic DNA elimination. The concentration and integrity of RNA were examined with a Thermo Scientific NanoDrop 8000 Spectrophotometer and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). RNA with OD260/280 ≥ 1.8 and RNA integrity number ≥7.0 was selected for the following experiment. Equal amounts of the high-quality RNA samples from gonad tissue were then pooled together for cDNA synthesis and sequencing.

Library construction, Illumina sequencing, and assembly.

RNA-Seq library preparation and sequencing was carried out by BerryGenomics Sequencing (Beijing, China). The cDNA library was prepared with ∼2.5 μg of total RNA based on the Illumina TrueSeq RNA Sample Preparation Kit (Illumina) protocols. The library was then amplified, and the final library yielded ∼400 ng with an average fragment size of ∼270 bp. The library was sequenced with one lane on an Illumina HiSeq 2000 instrument with 100 bp paired end reads after KAPA quantitation and dilution. Raw read data of yellow catfish RNA-Seq with SRA format have been uploaded to the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) under the accession number of SRR1103702. The clean reads from the four transcriptome filtered out adaptor-only reads and low-quality reads (reads with Q value ≤ 20). Cleaned reads were used for de novo assembly by the de Brujin graph approach with Trinity by default parameters.

Functional annotation and ontology.

The assembly RNA-Seq contigs were used for similarity search program against reference protein sequence including zebrafish, medaka, stickleback, and tetraodon. We used the BLASTX program with an e-value cutoff 1e-10 to perform similarity searches. Gene annotation was assigned to the RNA-Seq contigs based on the top BLAST hit. Gene ontology (GO) annotation analysis was then followed with Blast2GO. The level-two GO terms associated with transcriptome contigs were retrieved, and then the annotation result was categorized as biological process, molecular function, and cellular components. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were assigned to those assembled contigs using the online KEGG Automatic Annotation Server (http://www.genome.jp/tools/kaas).

Assembly assessment.

To compare the similarity to other teleost species, we compared the transcriptome contigs with RefSeq and Ensemble proteins of zebrafish, medaka, stickleback, and tetraodon, as well as the transcriptome of yellow catfish by using the BLAST program with default parameters.

Gene expression level calculation.

The read per kilobase of exon per million mapped reads (RPKM) values were calculated for the gene expression level with edgeR Bioconductor package release version 2.13. (http://bioconductor.org). The Kal's test was used to identify the differentially expressed genes with P value < 0.05.

Identification of specifically expressed genes, differentially expressed genes, and coexpressed genes.

Audic's method was used to identify differentially expressed genes (DEGs) between two libraries (1). The threshold for the P value was determined using false discovery rate (FDR) and was widely set at 0.05 (2). In this study, “FDR < 0.05” and “XX_RPKM = 0 or XY_RPKM = 0” were used to classified as specifically expressed genes (SEGs). “FDR < 0.05” and “|log₂(XX_RPKM/XY_RPKM)|> 1 or |log₂(XY_RPKM/XX_RPKM)| > 1” were used to classified DEGs. Those meeting “XX_RPKM < 2 and XY_RPKM < 2” statistical criteria were classified as nonexpressed genes (NEGs), whereas all the remaining ones were designated as coexpressed genes (COGs).This way, all genes were classified into four types: SEGs, DEGs, COGs, and NEGs.

Experimental validation by quantitative real-time PCR.

Selected differentially expressed genes was validated by quantitative real-time PCR (qRT-PCR) to verify the expression profile obtained from the transcriptome data. Gonads were dissected from XX and XY yellow catfish at one-year-old and two-year-old developmental stages. The total RNA was isolated from each sample and reverse-transcribed using MMLV reverse transcriptase (Invitrogen, Carlsbad, CA). The qRT-PCR Primers were designed using Primer Premier 6. The SYBR Green I Master Mix (TaKaRa, Dalian, China) was used for qRT-PCR using an Applied Biosystems Prism 7500-fast real-time PCR system. The PCR reactions were initiated by denaturation at 95°C for 5 min, followed by 40 amplification cycles at 95°C for 15 s and 60°C for 30 s. Dissociation protocols were used to measure melting curves and to control nonspecific signals from the primers. Results were expressed relative to the expression levels of β-actin in each sample using the Relative Expression Software Tool (REST) version 2009. The biological replicate fluorescence intensities of the control and treatment products for each gene, as measured by crossing-point values, were compared and converted to fold differences by the relative quantification method. Expression differences between groups were assessed for statistical significance using a randomization test in the REST software. The mRNA expression levels of all samples were normalized to the levels of β-actin gene in the same samples. Test amplifications were conducted with ensure that β-actin and target genes were within an acceptable range. A no-template control was run on all plates. qRT-PCR analysis was repeated in triplicate runs (technical replicates) to confirm expression patterns. The averages of three relative quantities of biological replications were used in a two-tailed Student's t-test with a 95% confidence level (P < 0.05) to determine the significance with respect to gene expression.

Transcriptome sequencing and assembly.

To obtain more detailed information about yellow catfish gonad transcriptome, four cDNA libraries derived from XX and XY gonads in 1-year juveniles and 2-year adults were sequenced by using Illumina HiSeq2000 technology. One-year and two-year old are noted as two developmental stage of yellow catfish in this study. As a result, a total of 211,701,078 paired-end reads of 100 bp were generated. The raw transcriptome sequences in this study have been uploaded in the NCBI SRA site; the accession number is SRR1103702. Low-quality sequences (Phred quality score < 30) and ambiguous nucleotides were removed. The remaining high-quality reads, 193,477,064 (91.4%), were obtained for transcriptome assembly and analysis (Table 1). The clean reads were assembled by using Trinity de novo assembling software (17). A total of 82,123 contigs were assembled, ranging from 351 to 21,268 bp, and the average length was 1611.2 bp, with an N50 value of 2,329 bp and an N90 value of 712 bp (Table 1). The transcriptome data analysis workflow is shown in Fig. 1.

Gene identification and functional annotation.

BLAST-based gene identification was performed to annotate the yellow catfish transcriptome and identify sex-related genes. BLASTX was used for the assembled contigs gene identification by searching against the reference protein sequences, including zebrafish (Danio rerio), medaka (Oryzias latipes), stickleback (Gasterosteus aculeatus), and tetraodon (Tetraodon nigroviridis). A total of 45,390 contigs had significant hit, corresponding to 21,869 unigenes. To further understand the orthologous gene clusters between yellow catfish and four other teleost fish species (zebrafish, medaka, stickleback, and tetraodon) with genome sequence information, a Venn diagram was drawn based on orthologous gene clusters among them (Fig. 2). We compared 19,588 yellow catfish gene families attained in this study with zebrafish (21,247 gene families), medaka (16,548 gene families), stickleback (17,354 gene families), and tetraodon (16,060 gene families); 7,988 gene families were shared among five teleost fish species (Fig. 2).

The differentially expressed unigenes were then used as input to perform GO annotation by Blast2GO (8). A total number of 26,804 assembled contigs were assigned at least one GO term, corresponding to 10,925 unigenes for describing biological processes, molecular functions, and cellular components (Table 1). The GO annotations were plotted in Fig. 3. Of these, the biological process ontology was the most prevalent, followed by the cellular component ontology and the molecular function. In brief, for biological processes, genes involved in cellular processes and single-organism processes were highly represented; for the cellular component, cell and cell part were the most represented categories; binding was the most represented GO term, followed by catalytic activity for molecular functions. Interestingly, within biological processes, a large fraction of unigenes are identified as related to growth, sex dimorphism, and reproduction. Among them, a total of 757 unigenes were annotated to reproduction (GO: 0000003), 5,390 unigenes were annotated to developmental process (GO: 0032502), 699 unigenes were annotated to growth (GO: 0040007), and 702 unigenes were annotated to reproductive process (GO: 0022414). In addition, the significantly enriched unigenes were annotated to regulation of developmental process (GO: 0050793) in 2-year-old adult testis comparing with ovary tissue in our study.

Contigs with previously identified gene matches were carried forward for further analysis. Function pathways analysis was based on the KEGG database, as an alternative approach for functional categorization and annotation for the differently expressed genes. Enzyme commission numbers were assigned to 18,600 sequences, which categorized them into different functional groups (Table 2). In brief, of the sequences with KEGG annotation, 4,691 (25.2%) were classified into metabolism, including major subgroups of global and overview maps (1,790 sequences), carbohydrate metabolism (566 sequences), lipid metabolism (471 sequences), and amino acid metabolism (455 sequences). Sequences grouped into genetic information processing accounted for 7.7% (1,435 sequences), including folding, sorting, and degradation (517 sequences); translation (468 sequences); replication and repair (250 sequences); and transcription (200 sequences). Environmental information processing, cellular processes, and organismal systems groups contained 3,497 (18.8%), 2,425 (13.0%), and 6,552 (35.2%) KEGG annotated sequences, respectively.

Assessment of transcriptome assembly.

To evaluate yellow catfish transcriptome assembly result, we compared the assembled contigs of yellow catfish transcriptome with RefSeq proteins of zebrafish, medaka, stickleback, and tetraodon by using the BLASTx program with an e-value cutoff of 1E-10. A total of 48,872 (59.5%), 41,041 (50.0%), 42,354 (51.6%), and 40,394 (49.2%) significant contig hits were identified with zebrafish, medaka, stickleback, and tetraodon, respectively (Table 3). After filtering the repeat isogenes, we found 26,239 unigenes in zebrafish, 24,661 in medaka, 27,576 in stickleback, and 23,118 in tetraodon (Table 3). In addition, 16,023 unique protein hits were identified, with zebrafish differing noticeably from those with medaka (12,917), stickleback (13,525), and tetraodon (12,830) (Table 3). Zebrafish is particularly striking in similarity to yellow catfish at the gene expression level. However, the similarity was still lower than what we expected. The transcriptome data presented here do not cover the whole genetic information due to limited tissue collection. To facilitate better understanding and characterization of yellow catfish transcriptome, more transcriptome data would be necessary from more tissue across all the crucial life stage and every circumstance or the whole-genome sequencing and assembly.

Microsatellites and structure variation identification.

For further assessment of the assembly quality and development of new molecular markers, we initially identified a total of 29,903 microsatellites from 82,123 contigs, including di-, tri-, tetra-, and pentanucleotide simple sequence repeats (SSRs) with a minimum of four repetitions for all motifs. The SSRs that were only located in one single read had been filtered. Among the 29,903 microsatellites, dinucleotide repeat motifs were the most abundant type (17,045, 57.0%), followed by trinucleotide repeat motifs (36.4%), tetra- (5.9%), and penta- (0.7%) (Table 4). These SSRs provide plenty molecular information to design polymorphic primers for further genotyping analysis.

For extended application of the RNA-Seq data, structure variations were discovered using the assembled transcriptome. The short reads of RNA-Seq data were aligned onto the reference transcriptome of yellow catfish; the total number of INDEL (insertion and deletion) varied from 4,687 to 8,163 from four different yellow catfish samples. The insertion number varied from 3,268 to 5,361, and the deletion number varied from 1,419 to 2,802. The total number of single nucleotide polymorphisms (SNPs) varied from 28,234 to 60,470 from four different yellow catfish samples. The synonymous number varied from 2,869 to 8,047, and the nonsynonymous number varied from 20,187 to 57,228 (Table 5). The available microsatellites and structure variants developed in this study serves as precious molecular markers for yellow catfish genetics research.

Sex-related gene expression profiling in gonads.

Of the 21,869 unigenes from yellow catfish transcriptome, 9,919 were found to be coexpressed in both XX and XY gonads, and 51 and 229 genes were detected to be expressed exclusively in XX and XY gonads in 1-year juveniles, respectively (as shown in Table 6 and Fig. 4). For 2-year-old adults, 7,911 genes were found to be coexpressed in both XX and XY gonads, and 40 and 1,188 genes were detected to be expressed exclusively in XX and XY gonads, respectively.

The number of XY-SEGs from two developmental stages was much greater than XX-SEGs. The number of SEGs was negatively correlated with the number of COGs expressed at two stages. While the numbers of XX-DEGs at ages 1 and 2 yr were 3,181 and 3,298, respectively, the corresponding numbers for XY-DEGs were 3,456 and 4,989, respectively (Table 6). Overall, more SEGs and DEGs were observed in XY gonads. A simple comparison of the scatter plots of the gene expression profiles at two developmental stages also revealed that there were more upregulated genes in XY than in XX (Fig. 4). Despite the significant differences in gene expression observed between XX and XY gonad transcriptome, a preliminary analysis revealed that the COGs expressed at two stages made up a significant proportion of the total (Table 6).

Identification of DEGs.

Identifying sex DEGs in 1 yr juveniles and 2 yr adults is regarded as an important approach in detecting molecular differences regulating sex determination and sexual dimorphism from the gonad transcriptome. Within two development stages, sex DEGs were extracted based on the FDR value (FDR < 5%) (Tables 7 and 8).

XY-DEGs in 1-year-old testis tissue (noted as YC1XY), were found, including FSHR (5, 30), DMRT1, LHR (45), CYP17a (52), LEPROT (16), GH, brain-derived P450 (20), EsR2, HSP70 (42, 55), Sox9a2, IGFBPII (66), and RPL15. Compared with YC1XY, more XY-DEGs were detected in the 2-year-old yellow catfish testis tissues (noted as YC2XY). Besides those genes previously identified in YC1XY, additional genes were detected particularly enriched expressed in 2-year-old yellow catfish, including GAPDH, CYP1a, 20BHSD, Sox9a1, GTHa. However, only two genes were expressed biasedly in ovary tissue during two developmental stages, which were GDF-9 and MyoD. Those genes with identities are known to be involved in gonadogenesis, spermatogenesis, testicular/ovarian development and differentiation, and sex determination genes could be identified among these sex-related genes.

Identification of SEGs.

The identification and characterization of the XX- and XY SEGs at two developmental stages is of vital importance for understanding the underlying molecular differences regulating XX and XY gonadal differentiation and development. Some genes were found exclusively involved in XY gonad development. For 1-year-old juveniles, only four sex-related genes uniquely expressed in XY gonads were identified: 1) rhophilin-associated protein 1 (ROPN1), an important reproduction-related gene, has been identified as a spermatogenic cell-specific protein and may be involved in sperm maturation, motility, capacitation, hyperactivation, and acrosome reaction (7, 13, 34); 2) sperm-associated antigen 6 (SPAG6), which has been shown to be a critical protein in either the assembly or structural integrity of the sperm tail axoneme (28, 68, 69); 3) sperm-associated antigen 16 (SPAG16), which plays an essential roles for spermatogenesis, germ cell viability, and the integrity of the axoneme (67); and 4) fibroblast growth factor-binding protein 2 (FGFBP2); it has been found that high FGFBP2 expression in high-grade gliomas is positively correlated with survival, and it also is more closely correlated with survival than histological grade (11, 62) (Table 9).

We also identified the following set of 10 genes exclusively expressed in 2-year-old adults: 1) SPAG6; 2) TSSK [the members of the TSSK family may have a role in sperm differentiation in the testis and/or fertilization. Also, TSSK have an important role in germ cell differentiation and/or sperm function (37, 47, 50, 51)]; 3) dynein heavy chain, including 1, 2, 3, and 12 (23); 4) ANKMY1; 5) SPAG16; 6) TCTE1, which is absolutely required for fertilization and expressed in earlier stages of spermatogenesis and also may facilitate species-specific divergence of sperm function (36); 7) dual-specificity testis-specific protein; 8) TESK (several studies revealed that TESK mRNA is highly expressed in testicular germ cells at the stages of pachytene spermatocytes to round spermatids); 9) sterol 26-hydroxylase (4); and 10) androgen receptor (AR) [androgens are necessary for spermatogenesis in most studied vertebrates, including teleosts (22)]. These genes were only expressed in male yellow catfish (Table 10); their functions are certainly worthwhile of additional investigative effort. As male-specific genes in yellow catfish, they represent interesting candidate transcripts for experimental validation by PCR amplification.

Validation of RNA-Seq results by qPCR.

To validate the DEGs identified by RNA-Seq, we randomly selected nine genes from those with differential expression patterns and from DEGs and SEGs based on their function for qRT-PCR validation. Comparing the transcriptome data with the qRT-PCR results from nine selected DEGs, we found the qRT-PCR results to be consistent with the RNA-Seq results at two developmental stages (Fig. 5). Overall, the qRT-PCR results indicate the reliability and accuracy of the Trinity reference assembly and the RNA-Seq-based transcriptome expression analysis.

FSHR/LHR.

We identified FSHR/LHR as XY DEGs. As glycoprotein hormones, follicle-stimulating hormone (FSH) and luteinizing hormone (LH), called gonadotropins, are believed to regulate gonadal functions (35). Gonadotropins exert their action through gonadotropin receptors, the LH receptor (LHCGRBB) and the FSH receptor (FSHRA). Because of their importance in the regulation of reproduction, including ovulation, they have been studied in many species. FSH/FSHRa regulates the early phases of gametogenesis, such as vitellogenesis, whereas LH/LHCGRbb stimulates the final maturation stages, such as ovulation (33, 35, 44, 64). In our findings, the fact that both FSHR and LHR were identified among DEGs in XY gonads is intriguing, suggesting that they may have similarly functions in regulating reproduction.

Dmrt.

Of great interest in our study was detection of highly expressed Dmrt1 as an XY-DEG belonging to the Dmrt gene family. Dmrt genes have been reported to be actively involved in sex determination and/or differentiation. These genes stimulate male-specific differentiation but repress female-specific differentiation (21, 29, 60). Dmrt1 is claimed as a master sex-determination gene in medaka (41). In a report on Siberian sturgeon, Dmrt1 showed a significantly higher expression in testis. A similar trend of dmrt1 expression during early and advanced stages of gonad development has been observed in sturgeon (3), rainbow trout (40), catfish (53), and zebrafish (18). In this study, there was significant Dmrt1 gene expression in testis tissues in both developmental stages.

Sox9.

The abundance of expression of Sox9a2 and Sox9a1, as members of Sox 9, was examined in yellow catfish as XY-DEGs. Sox9 is a transcriptional factor that is a member of the Sox gene family of SRY-related HMG box genes (12). Previous studies have indicated that Sox9 plays a very important role in the male gonadal development of many vertebrates (6, 27, 43, 46, 56, 58, 59). Sox9 as a critical factor in sex determination and differentiation was also confirmed by recent studies in teleost fish. It is found expressed in testis of rainbow trout (54) and in channel catfish (53). Sox9 expressed during testicular developmental stage in fish is considered as a candidate gene involved in testis differentiation but not in sex determination.

Other genes.

Insulin-like growth factor/insulin-like growth factor-binding protein (IGF/IGFBP), exhibiting high levels of expression in yellow catfish gonads, were identified as DEGs in XY gonads. The IGF system plays a pivotal role in vertebrate growth, development, proliferation, and metabolic regulation. Some genes are believed to be critical to gonadal development in both sexes, such as IGFBP (65) and IGFBP2 (66). IGFBP2 has been posited to have important functions in growth, development, and other physical actions in Japanese flounder (66). Other studies also indicate that igfbp2 could regulate the growth of zebrafish by acting downstream in the growth hormone (GH)/IGF-I axis (7, 10, 14, 25).

LEPROT was found to be significantly highly expressed in testis compared with the ovary in yellow catfish, and its gene expression pattern is consistent with findings by Gong et al. (16). Highly expressed LEPROT was also detected in testis of other fish species (31, 32, 39, 57, 61). In teleost fish, cytochrome P450 aromatase (P450arom) is predominantly expressed in the central nervous system and in the gonads. In our study, P450 was detected as XY-DEGs in gonads.

As for XX-DEGs, growth and differentiation factor 9 (GDF9) and MyoD were highly expressed in ovary in our study. GDF9 is a key regulator of follicular development and might have a determining role in establishing the ovulation quota, which has raised an entirely new set of questions regarding the control of follicle/oocyte maturation (15, 24, 48).

In vertebrates, MyoD is one of the important myogenic regulatory factors (MRFs) that are directly or indirectly regulated by growth hormone. Evidence indicates that muscle growth is regulated positively and negatively by GH directly or indirectly act on MRFs. Detailed examination for those genes should be directed toward understanding their function in sex determination of yellow catfish.