Comparative Transcriptome Analysis of Differentially Expressed Genes and Signaling Pathways between XY and YY Testis in Yellow Catfish

YY super-males have rarely been detected in nature and only been artificially created in some fish species including tilapia and yellow catfish (Pelteobagrusfulvidraco), which provides a promising model for testis development and spermatogenesis. In our previous study, significant differences in morphology and miRNA expression were detected between XY and YY testis of yellow catfish. Here, solexa sequencing technology was further performed to compare mRNA expression between XY and YY testis. Compared with unigenes expressed in XY testis, 1146 and 1235 unigenes have significantly higher and lower expression in YY testis, respectively. 605 differentially expressed unigenes were annotated to 1604 GO terms with 319 and 286 genes having relative higher expression in XY and YY testis. KEGG analysis suggested different levels of PI3K-AKT and G protein-coupled receptor (GPCR) signaling pathways between XY and YY testis. Down-regulation of miR-141/429 in YY testis was speculated to promote testis development and maturation, and several factors in PI3K-AKT and GPCR signaling pathways were found as predicted targets of miR-141/429, several of which were confirmed by dual-luciferase reporter assays. Our study provides a comparative transcriptome analysis between XY and YY testis, and reveals interactions between miRNAs and their target genes that are possibly involved in regulating testis development and spermatogenesis.

As an evolutionary link between invertebrates and higher vertebrates, fish species have a very complex sex determination system with XX/XY male heterogametic system as the main sex determination system [16]. Until now, YY super-males have only been artificially produced in some fish species [17][18][19], but they could not survive in case of higher vertebrates including human and mouse. Previous studies on rainbow trout (Oncorhynchusmykiss) reported obvious differences in gene expression and morphology between XY and YY testis [20,21]. Significant differences of aromatase expression were found in spermatogonia, spermatids and epithelial cells among XY and YY testis [21]. Moreover, a relative higher level of androgen receptor expression was observed in efferent ducts of YY testis compared with XY testis [20]. However, the sperm quality and quantity of XY and YY males were the similar in Nile tilapia (Oreochromisniloticus) [22]. There are no comprehensive studies regarding genetic differences between XY and YY males.
In our previous study, we observed larger spermatogenic cyst and more spermatids in YY super-males than XY males by histological analysis, suggesting a higher degree of testis maturity. Many miRNAs that are potentially involved in testis development and spermatogenesis were identified to be differentially expressed between XY and YY testis [23]. The establishment and maintenance of spermatogenesis in fish requires specialized gene regulatory networks in the testis [24]. Here, we utilized RNA-Seq approach to identify genes and pathways that were differentially expressed between XY and YY testis and their functional relationship with miR-NAs. Hopefully, our findings would provide a clue about the genetic mechanism of testicular differences between XY and YY males.

Illumina sequencing, sequence assembly and functional annotation
In order to identify differentially expressed mRNA in testes of male and super-male yellow catfish, two solexa libraries were constructed by XY and YY testes respectively. After removing adaptors and low quality reads, a total of 87,149,414 and 75,448,188 reads were obtained in each profile respectively. After de novo assembly, 78148 unigenes were obtained by paired-end method of Trinity and TGICL clustering with mean length of 944bps (S1 Fig). The unigenes of de novo assembly were searched against the NCBI non-redundant (nr), SWISS-PROT, KEGG, GO and KOG protein databases by using BLASTX with a cut-off Evalue of 1e -5 (S1 Table). Finally, a total of 19,795 (25.33%) unigenes were significantly matched with nr database. Among these BlastX-hit unigenes (Fig 1), 11591 (58.56%) hits were assigned to Daniorerio while only 776 (3.92%) hits overlapped with Ictalurus punctatus (the close species to Pelteobagrus fulvidraco), perhaps due to the limited amount of genomic data on Gene-Bank about the species of Siluriformes. In addition, there were a number of unigenes without hitting results, which may contain non-coding fragments or some undetected genes.

Analysis of differentially expressed genes (DEGs)
To identity differentially expressed unigenes between XY male and YY super-male, the expression of assembled unigenes were counted by RPKM method. After comparison (XYvsYY testis) with fold change threshold value = 2 and FDR test (P< 0.05), 4458 unigenes were found expressed with significant difference. Among them 1006 were detected from testis of YY supermale and 1072 were found from testis of XY male. Besides this, 1146 unigenes from YY and 1235 from XY testis were detected with significantly higher expression. The MA scatter plots comparison reveal the differential expressed genes existed between testes of male and supermale in yellow catfish (Fig 2).

GO annotation and analysis of DEGs enriched in YY testis
The 4458 differentially expressed unigenes (DEGs) were searched against the Gene Ontology database (www.geneontology.org) to determine which kinds of GO term DEGs mainly participated in. Finally, 605 DEGs were annotated to 1604 GO terms (S2 Table) accompanied by 319 and 286 genes with relative higher expression in XY and YY testes, respectively. The functions of 286 DEGs with higher expression in YY than XY testis were assigned to biological process, cellular component and molecular function (Fig 3). In biological process, proteolysis (58 DEGs, GO:0006508), RNA-dependent DNA replication (38 DEGs, GO: 0006278) and DNA integration (33DEGs, GO:0015074) were the most prominent terms. Integral component of membrane (112 DEGs,GO:0016021) was the most prominent within cellular component followed by nucleus (67DEGs,GO:0005634), cytoplasm (48DEGs,GO:0005737) and plasma membrane (40DEGs, GO:0005886). In molecular function, most of the annotated unique sequences were assigned to zinc ion binding (74DEGs,GO:0008270), ATP binding (51DEGs,GO: 0005524), metal ion binding (53DEGs,GO:0046872), RNA binding (46 sequences,GO:0003723) and RNA-directed DNA polymerase activity (38DEGs,GO:0003964) (Fig 3).
The PI3K-AKT signaling pathway is involved in many fundamental functions including testis development and spermatogenesis, and stimulated by many kind of growth factors that specifically binds to receptor tyrosine kinase (RTK) or G protein-coupled receptors (GPCR) [25,26]. As members of RTK, the spleen tyrosine kinase (Syk), colony stimulating factor 1 receptor (Csf1r), and prolactin receptor (Prlr) were expressed about 3.68, 3.71 and 1.73 fold higher in YY testis than in XY, respectively. In addition, β1-Integrin (Itgb1) and β2-Integrin (Itgb2) were also expressed 2.71 and 3.15-fold higher in YY testis than in XY (Fig 5). In the "neuroactive ligand-receptor interaction" signaling pathway, expression of multiple genes associated with G protein signaling were significantly up-regulated in YY testis (Fig 6), like Kiss1r (GPR54) and somatostatin receptors (Sstr) in Class A of GPCR signaling and metabotropic glutamate receptor 5 (GRM5) in Class C of GPCR signaling. Meanwhile, the mRNA levels of glutamate receptor AMPA 2b (gria2b), glutamate receptor AMPA 4a (gria4a) and prolactin receptor (PRLR) were also higher in YY than in XY. In contrast, the expression levels of histamine receptor H1 (Hrh1), neuropeptide FF receptor 2 (Npffr2), leukotriene B4 receptor 1 (Ltb4r1) and nicotinic acetylcholine receptor α1 (Chrnα1) were higher in XY than in YY testis.

qRT-PCR confirmation of DEGs between XY and YY testis
To verify the accuracy of the sequencing data, twelve DEGs related to RTK and G protein signaling pathway were arbitrarily selected and validated by quantitative real-Time PCR (qRT-PCR). The twelve DEGs includes four genes (Hrh1, Npffr2, Ltb4r1, Chrnα1) relatively high expressed in XY and eight genes (Prlr, Csf1r, Itgb1, Kiss1r, Sstr, GRM5, Gria2b, Gria4a) relatively high expressed in YY testis in the transcriptome data ( Fig 7A). As in the qRT-PCR results, the relative expression levels of twelve differentially expressed genes were completely consistent with the Solexa sequencing ( Fig 7B).
Identification of potential target genes for miR-141/429 from DEGs and validation by dual-luciferase reporter assays In our previous study, relative lower expression of miR-141/429 was observed in YY testis indicating a higher degree of testis maturity than XY testis. A high dose of 17α-ethinylestradiol (EE2) up-regulates the expression of miR-141/429 [23]. Here, we examined whether some DEGs are potential target genes for miR-141-3p and miR-429b-3p. Finally, 31 and 11 YY enriched DEGs were predicted to be targeted by miR-141-3p and miR-429b-3p (Table 1). For example, Itgb2, a factor involved in the PI3K-AKT signaling pathway, was highly expressed in YY and was a putative target of miR-141-3p. In addition, Gria4a is a factor for the neuroactive ligand-receptor interaction signaling pathway and also a putative target for miR-141-3p. AMH and Tgfβr1 were potential targets for miR-141-3p and miR-429b-3p, respectively.
To determine whether there are direct interactions between miR-141-3p and PI3K-AKT or GPCR signaling pathway, we used dual-luciferase reporter assays to measure the inhibitory effect of this miRNA on Itgb2 and Gria4a. There is one binding site for miR-141-3p in either Itgb2 or Gria4a 3'UTR ( Fig 8A). We sub-cloned the 3'UTR of Itgb2 or Gria4a into the pmir-GLO vector, and co-transfected each construct with miR-141 mimic or its appropriate control into HEK293 cells. The results showed that miR-141-3p down-regulated luciferase activity by 69% (±3%) in Gria4a 3'UTR and 26% (±4%) in Itgb2 3'UTR, respectively (Fig 8B). These results strongly support the prediction of Itgb2 or Gria4a as direct targets of miR-141-3p.

Discussion
In vertebrates with XY sex-determining system, the expression of sex-determining gene on Y chromosome leads to the development of male phenotypes and testis. Yellow catfish, an important aquaculture fish in China has XY sex-determining system and displays sexual size dimorphism as male grows 2-3 times faster than female. 454 pyrosequencing and illumina sequencing studies have been performed to compare differentially expressed genes and pathways between XX ovary and XY testis and provide a valuable genomic resource for studying fish reproduction, sex determination and differentiation [27,28]. However, there are limited studies regarding the gene expression difference between XY and YY fish. Therefore, we used solexa sequencing technology to compare mRNA expression between XY and YY testis of yellow catfish.
In fully mature 18 month-old Nile tilapia, the sperm quality and quantity of XY and YY males were the similar [22]. However, Herrera et al. found that YY tilapia has superior reproductive capacity than XY fish, since the primordial germ cells and spermatogenic cells in YY were larger than XY fish during gonad development and the lobules, blastemal of the reproductive duct and mature sperm cells appeared earlier in YY than XY fish. Moreover, YY tilapia has bigger testis, thicker somatic tissues and more spermatogenic cells, as well as matures earlier than XY fish [29]. In yellow catfish, we also found that YY matures earlier and has superior reproductive capacity than XY fish, as larger spermatogenic cyst and more spermatids were observed in YY fish [23]. Interestingly, androgen receptor has a relative higher level of expression in efferent ducts of YY testis compared to XY testis in rainbow trout [20]. All these studies suggest that YY testis has substantial difference in histology, structure and gene expression compare to XY testis.
It is important to know the context under which specific signaling pathways regulate sperm maturation as well as testis development in YY that matures earlier than XY fish. The PI3K-AKT signaling pathway plays essential roles in testis development and spermatogenesis, as loss of p110beta subunit of phosphoinositide 3-OH kinase impaired spermatogenesis and lead to defective fertility [26]. Activation of the PI3k/Akt pathway by membrane progestin receptor-alpha stimulated sperms leads to hypermotility in Atlantic croaker [25]. Lgr4, one of the orphan GPCRs regulates sperm development and fertility [30]. Testosterone signaling is mediated by a G-protein-coupled receptor and its interactors [31]. We have found more PI3K-AKT and GPCR signals in YY than XY yellow catfish, such as syk, prlr and kiss1r, coincides with a higher degree for testis maturity and advantageous spermatogenesis in YY than XY yellow catfish [23]. As one of the cytoskeletal component of spermatic flagella, tyrosine-  phosphorylated Syk could bind and phosphorylate to its downstream part PLCγ1 and regulate the metabolic hyperactivated motilityof spermatozoa [32][33][34].
In mammals, increased expression of miR-141/429 was associated with defects of spermatogenesis [35][36][37]. High dose of 17α-ethinylestradiol (EE2) resulted in upregulation of miR-141/ 429 and impairment of spermatogenesis [23]. In our study, miR-141/429 was predicted to target several factors in PI3K-AKT and GPCR signaling pathways that were involved in testis development and spermatogenesis. Further characterization of the interaction of miRNAs and their targets could contribute to a better understanding about the molecular mechanisms of testis develop and spermatogenesis.

Solexa library construction and sequencing
All experimental procedures involved in this study were permitted by the Institutional Animal Care and Institute of Huazhong Agricultural University. The total RNAs of XY (4 individuals) and YY (3 individuals) testes were the same biological sample as described before [23]. The RNA integrity analysis was performed by the Agilent 2100 Bioanalyzer (Agilent Technologies). An amount of 3μg total RNA was used for libraries construction. All the following process was performed by using the Illumina RNA Sample Preparation Kit, following the manufacturer's protocols. The mRNA was concentrated by oligo(dt) magnetic beads and then made into short fragments (~200bp) to work as templates for synthesizing the first strand cDNA. The double strand cDNA libraries were synthesized and purified by agarose gel purification, in which, the cDNA fragments were coupled by sequencing adptorsat the 5' and 3' ends. The male (XY) and super-mal (YY) libraries were sequenced on an IlluminaHiSeq 2000 platform. After removing the adaptor, low quality bases (Length threshold value>35bp), 3'-end low quality bases (Quality threshold value>20) and the reads containing the "N", the clean reads were assembled into contigs with the Trinity software [30]. The generated contigs were clustered into unigenes by performing TGICL software system [31]. The raw reads of yellow catfish gonad have been deposited to the NCBI database (accession no: SRR1845493).

Quantification of differential expressed unigenes
To detect the differential expression of unigenes, DESeq software package with RPKM (Reads Per Kb Million reads) method was performed to quantify the expression of two expression profiles. The formula applied was RPKM ¼ 10 6 C NL=10 3 , in which C is the number of reads mapper to merged transcripts, N and L are the total mapped reads and base number of one unigene [38]. The False Discovery Rate (FDR) method was applied in multiple hypothesis of test to correct significant levels and eliminate influence of random fluctuations and errors [39]. After calibration, the ratios of RPKMs were used to calculate fold-change with threshold value cut-off of 2-fold and Negative binomial distribution hypothesis-testing with P< 0.05 [40].

Functional annotation and GO/KEGG enrichment analysis
The unigenes were searched against databases of NCBI nr, SWISS-PROT, TrEMBL, Cdd, pfam and KOG by Blast X with a cut-off E-value of 10e -5 . The results of BlastX annotation were uploaded on Blast 2 Go to generate Gene Ontology annotations and mapped to the categories of GO database ((geneontology.org/page/download-annotations) [41]. The results of BlastX were also searched against the Kyto Encyclopedia Genes and Genomes (KEGG) in Blast 2Go. All the differential expressed unigenes were mapped to KEGG database with counted numbers of involved interactive metabolic pathways (p<0.05). To investigate which GO item and signaling pathways the DEGs participated in, all of the clustered DEGs were mapped back to the GO and KEGG databases. Each enrichment item was corresponded to a specific enrichment score that was calculated by the formula, score value = m n =M N . The statistical significance of the GO enrichment was evaluated by the hypergeometric distribution testing, P = 1- where N is the number of unigenes with GO annotation, n is the number of DEGs with GO annotation, M is the number of unigenes with one specific GO annotation and m is the number of differently expressed unigene with one specific GO annotation [42].

Quantitative real-time PCR (qRT-PCR)
Briefly, 1μg of total RNA was reverse transcribed by using PrimeScriptRT reagent Kit (Takara) according to the protocol. The qRT-PCR reaction was performed in a 20μl reaction volume using the Roche LightCycler 480 with SYBR Green PCR master mix(Roche) and three biological replicates were conducted for each reaction. The β-actin was selected as an internal reference to normalize the Ct values of each reaction by using the 2 -ΔΔCt method [43]. The ANOVA analysis was used to perform differential expression analysis.

Identification of direct miRNA targets and validation by dual-luciferase reporter assays
In our recent study, we found that miR-141-3p and miR-429b-3p have higher expression level in XY than YY testis [23]. To investigate potential targets of miR-141-3p and miR-429b-3p, first the Open Reading Frame (ORF) and 3'UTR of YY highly expressed unigenes were predicted, by searching against the vertebrate genomic database in GENSCAN (http://genes.mit. edu/GENSCAN.html). Perl scripts of both Target Scan and miRanda were performed for searching the putative targets with default parameters, including context score percentile !100 for Target Scanand Max_Energy −20 and for miRandabased on hybrid energy and stability [44,45] To characterize the interaction between miR-141-3p/-429b-3p and their predicted target genes, the 3'-UTR of selected putative target genes (Itgb2 and Gria4a) were inserted into the pmirGLO expression vector (Promega, USA). Hek-293T cells were seeded in 96-well plates and co-transfected with the constructed vectors and microRNA mimics or its control oligonucleotide using DharmaFECT transfection reagent (Dharmacon). 36h post transfection, the dual-luciferase reporter assay system was used to detect reporter (Firefly and Renilla) activity as described [46]. The profile of relative luciferase activities were normalized to Renilla luciferase activities.