Gonadal Transcriptome Analysis of Male and Female Olive Flounder (Paralichthys olivaceus)

Olive flounder (Paralichthys olivaceus) is an important commercially cultured marine flatfish in China, Korea, and Japan, of which female grows faster than male. In order to explore the molecular mechanism of flounder sex determination and development, we used RNA-seq technology to investigate transcriptomes of flounder gonads. This produced 22,253,217 and 19,777,841 qualified reads from ovary and testes, which were jointly assembled into 97,233 contigs. Among them, 23,223 contigs were mapped to known genes, of which 2,193 were predicted to be differentially expressed in ovary and 887 in testes. According to annotation information, several sex-related biological pathways including ovarian steroidogenesis and estrogen signaling pathways were firstly found in flounder. The dimorphic expression of overall sex-related genes provides further insights into sex determination and gonadal development. Our study also provides an archive for further studies of molecular mechanism of fish sex determination.


Introduction
According to Helfman et al., there are almost 30,000 species of fish distributed nearly in all the aquatic habitats around the world [1]. They are the most abundant vertebrates on Earth, showing a diversity of species unmatched by other classes. Not surprisingly given this extreme diversity, fish exhibit all known forms of vertebrate sex determination to adapt to the variable habitats [2]. In the meantime, economic values of growth rate, time and age of maturation, body shape, and carcass composition are related to their sexual development in some edible fish species [3]. Then, we are interested in fish sex determination mechanism. As reported, fish sex determination patterns can be classified as genetic sex determination (GSD) and environmental sex determination (ESD) forms. However, the feature of sex determination in fish is remarkably flexible; even individuals with GSD can be influenced by environmental factors like temperature, that is, GSD + EE (environmental effects) form [4]. Although fish have several different sex determination forms, it is hypothesized that genes involved in sex determination are probably conserved throughout evolution. Several genes have been confirmed as master genes of sex determination in some fish species. In the medaka (Oryzias latipes), a homologue of the dmrt1 gene (called dmy) is located on the Y chromosome, and its expression is a necessary and sufficient condition for triggering testicular development in bipotential gonads [5,6]. Recently, five novel sex determining genes (or candidates) have been reported in other fish: amhy in Patagonian pejerrey (Odontesthes hatcheri) [7], irf9y in rainbow trout (Oncorhynchus mykiss) [8], gsdf in Oryzias luzonensis (a relative of medaka) [9], amhr2 in fugu (Takifugu rubripes) [10], and dmrt1 in half-smooth tongue sole (Cynoglossus semilaevis) [11]. Besides the sex determining genes, some conserved genes shown to play important roles in mammal sex determination and differentiation were cloned and identified in fish. These include cyp19, foxl2, sf1, dax1, wt1, mis, dmrt, and sox9 [12]. In mammals, these genes act together to constitute complicated network whereby sex phenotype is established [13]. However, studies on the function and connections of the above genes in fish are limited. More or novel sex-related genes are also needed to be found out. And then, the complex mechanism of fish sex determination could be adequately explained.
Over the past decade, significant progress has been made in genome-wide gene expression profiling by the development and application of large scale sequencing technique, which can easily show more differential expressional genes in different traits, such as gender. Transcriptome profiling associated with sex determination and differentiation using RNA-seq of several fish, including platyfish (Xiphophorus maculatus) [14], rainbow trout (O. mykiss) [15], Nile tilapia (Oreochromis niloticus) [16], rockfish (Sebastiscus marmoratus) [17], catfish (Ictalurus punctatus) [18], and turbot (Scophthalmus maximus) [19], was shown. These data provided transcriptomic information expressed in gonads at particular condition and time and identified sex differentially expressed genes, while there is almost no report on biological pathways including gonadal steroidogenesis pathway in male and female fish.
As an important commercially marine flatfish, olive flounder (Paralichthys olivaceus) is mainly cultured in China, Japan, and Korea. The female flounder grows significantly faster and bigger than the male one [20,21]. Breeding and culturing all-female population of flounder is a promising approach to boost production. Thus, the study on sex manipulation of flounder has been attracting researchers' interests. Olive flounder has XX (female)/XY (male) sex determination system indirectly inferred from the female-dominant phenotype among the gynogenesis offspring [21]. Its formation of sex phenotype is also influenced by environmental factors such as water temperature and external hormone [22]. Sexrelated conserved genes have been studied to elucidate the sexual molecular mechanism in flounder. Kitano et al. firstly cloned flounder cyp19a gene and found female predominant expression pattern [23], whereafter cloning and expression profile analysis of cyp19a and its transcription factors such as foxl2 [24] and dmrt1 [25] and other sex-related genes including cyp17 [26] and dmrt4 [27] were conducted with male and female flounder. However, we have got only limited message about flounder sex determination, and more sexrelated genes and their functions need to be studied. Further identification of the expression profile of genes involved in gonadal development using RNA-seq may help to illuminate the gene regulatory network controlling sex determination and subsequent maintenance of phenotypic sex. In this study, we sequenced flounder gonadal transcriptomes and identified the differences in gene expression profiles between ovary and testis and relevant biological pathways. These data would provide a useful genomic resource for future study on sex determination and for selection of candidate genes involved in these processes in flounder.

Fish.
Adult and juvenile flounder (12∼40 cm in total length, TL) used in the present study were collected from Shenghang fish farm (Weihai, China) or purchased from Nanshan market (Qingdao, China) and were temporarily cultured in a 3 m 3 aerated seawater tank at the institute aquarium and fed with commercial particle food twice a day. Gonads were retrieved from the abdomen of fish after anesthetization. Their genders were identified by morphological observation of gonads [28]. Each gonadal sample was divided into two halves. The first half was fixed in Davison's fixative solution for identification of gonadal developmental stage using a histology method as described below. The second half was stored immediately in liquid nitrogen for RNA isolation. Totally, the gonadal tissues at developmental stage I (3 testes and 3 ovaries), II (2 testes and 3 ovaries), and III (2 testes and 3 ovaries) were used in this project. Semen and eggs were gently squeezed out from 2 mature male flounder (40∼ 50 cm in TL) and 3 mature female flounder (50∼60 cm in TL), respectively, and the samples were immediately frozen and stored in liquid nitrogen for RNA isolation. All animal work has been conducted according to relevant national and international guidelines. Animal protocols were approved by the Institute of Oceanology, Chinese Academy of Science.

Gonadal Histology.
For histological analysis, ovary or testes tissue from each adult or juvenile flounder was fixed in Davison fixative solution for 24 h and stored in 70% ethanol. Pieces of gonadal tissue were cut down, dehydrated using ethanol (gradient: 70%∼100%), and finally embedded in paraffin. The sections were conducted at the thickness of 5∼7 m, dewaxed by ethanol (gradient: 100%∼50%), washed by distilled water, stained by hematoxylin-eosin, and observed microscopically after air drying. Slices of each sample were mounted on glass slides, stained with hematoxylin, and counterstained with eosin (HE staining) to determine the developmental stage of gonads. Detailed methods can be found in Sun et al. [22] and Radonic and Macchi [29].

RNA Isolation and cDNA Library Construction.
Total RNA was isolated from each sample using Trizol Reagent (Invitrogen, USA, http://www.lifetechnologies.com/ cn/zh/home/brands/invitrogen) based on the protocol. The genomic DNA was eliminated by treatment with DNase1 (10 U/mL, Ambion, USA, http://www.lifetechnologies.com/ cn/zh/home/brands/ambion) at 37 ∘ C for 1 h. The purified mRNA was enriched by Micropoly(A) Purist RNA purification kit (Amion, USA), and the concentration and integrity of mRNA were qualified using Agilent 2100 Bioanalyzer (Agilent Technologies, USA, http://www.home.agilent.com). The mRNAs from ovarian developmental stages I, II, and III and egg were mixed together to synthesize cDNA, and mRNAs from testicular developmental stages I, II, and III and sperm were mixed together correspondingly. These blended mRNAs served as templates to synthesize firststrand cDNA using GsuI-oligo dT primer; the reaction was performed with Superscript II reverse transcriptase (Invitrogen, USA) at 42 ∘ C for 1 h. Biotins were subsequently attached to the 5 cap of mRNA oxidized by NaIO 4 (Sigma, USA, http://www.sigmaaldrich.com), whereby the biotin-labeled mRNA /cDNA could be sublimed by Dynal M280 magnetic beads (Invitrogen, USA). Released from the hybrid strands by alkaline lysis, the first-strand cDNA was attached with an adaptor at its 5 end. Second-strand cDNA was synthesized using Ex Taq polymerase (Takara, Japan, http://www.takara-bio.com/) based on the first-strand modified cDNA; then poly(A) and 5 adaptor were trimmed by GsuI enzyme. The synthesized cDNA was disrupted into short fragments (300-500 nt) by ultrasound instrument which were further enriched by Ampure beads (Agencourt, USA, https://www.beckmancoulter.com). These purified cDNA fragments were used to construct cDNA library with the method of TruSeqTM DNA sample Prep kit-set (Illumina, USA, http://www.illumina.com/). Finally, the two cDNA libraries were sequenced on Illumina Solexa using paired-end strategy in a single run.

Illumina Sequencing, Functional Annotation, and Bioinformatics Analysis.
Total reads were produced through Illumina Solexa instrument (2 * 100 bp pair-end sequencing) from Chinese National Human Genome Center in Shanghai. The clean reads were obtained from original data by filtering out reads inclusive of unknown nucleotides and low-quality reads in which Q5 percentage (Q5 percentage is proportion of nucleotides with quality value larger than 5) is less than 50%. All clean reads of the two libraries were jointly assembled into contigs performed by Trinity software. The assembled contigs were conducted to predict protein-coding region by GetORF module of EMBOSS package [30]. All the protein-coding sequences were submitted for blastp similarity searches against the NCBI nonredundant (NR) protein database and Eukaryotic Ortholog Groups (KOG) database with the evalue of top hit lower than 1 e −5 . Furthermore, GoPipe software was used to perform blastp (cut-off e-value of <1 e −5 ) search against the Swiss-Prot database and TrEMBL database. With the result of blastp, gene ontology (GO) annotation associated with "biological process, " "molecular function, " and "cellular component" was obtained using the gene2go. Likewise, the predicted protein sequence was submitted for bidirectional blastp (cut-off e-value of <1 e −3 ) similarity searches against Kyoto Encyclopedia of Genes and Genomes (KEGG) database to assign KEGG Orthology (KO) number. According to the KO assignment, metabolic pathways were generated with tools supplied by KEGG [31]. The mapped read count of given gene is affected by its length and sequencing depth; the reads per kb per million reads (RPKM) were calculated to standardize gene expression level [32]: Here, indicates the mapped read count of a given gene from a given library. indicates the length of a given gene. indicates total mapped read count of a given library.

Identification of Sex-Related Differentially Expressed
Genes. RPKM was directly used to compare the difference of gene expression level between male and female. This process was completed by DEGseq (an package) based on the MARS model (MA-plot-based with random sampling model) [33]. We used Benjamini-Hochberg method to determine the threshold of the -value in multiple testing. In our study, "q < 10 −3 " and "|log 2 (RPKM XX/RPKM XY)| ≥ 2" were chosen to identify sex-biased genes. Furthermore, we adopted four strategies to excavate sex-related genes from them: (i) complete cDNA sequences of well-known sex-related genes were downloaded from NCBI nucleotide database (http://www.ncbi.nlm.nih.gov/nucleotide/) and were conducted local blast (cut-off e-value of <1 e −10 ) search against the local contigs database; (ii) sex-related keywords were retrieved in the annotation of sex-biased genes; (iii) among the sex-related KEGG pathways, the particular genes encoding pivotal enzyme were selected; (iv) genes reported to be relevant to sex differentiation were selected from the sexbiased genes.  (Table S1). The function of the predicted proteins was classified with GO assignments statistically analyzed in Figure  S3. In total, there were 20,582 proteins assigned with 132,920 GO terms and the three corresponding organizing principles.  (Table S1 KEGG pathway). Figure S4 showed the category distribution of biological pathways, in which the most enriched pathway is associated with signal transduction.

Sex-Biased Genes.
Total 10348 male-biased and 3296 female-biased contigs were identified, respectively, (Table  S2) showing significant expression difference between male and female. Among these sex-biased contigs, 887 (8.57%) male-biased contigs and 2193 (66.54%) female-biased contigs were annotated with known genes. In addition, 2111 contigs were identified to be expressed specifically in male and 75 contigs were identified to be expressed specifically in female. The statistical overview of transcriptome analysis was shown in Table 1. By using four strategies, sex-related well-documented genes were identified. In strategy i, we searched for well-known candidate genes already characterized to be sexually dimorphic in flounder. Among these genes, SRY-box containing protein 9 gene (sox9), sox8a mullerian inhibiting substance (mis), doublesex and mab-3 related transcription factor 1 gene (dmrt1), and so forth were predominantly expressed in male library, whereas P450 aromatase gene (cyp19a), forkhead transcription factor L2 gene (foxl2), orphan nuclear receptor dax1, and so forth were obviously overexpressed in the female library. We also found out sox6b gene containing SRY-box, which is not reported in fish.In strategy ii, a set of keywords, including male, female, sex, sperm, egg, ovary, testis, estrogen, and androgen, were used to search sex-related genes based on annotation results. Amounts of genes consisting of the above keywords such as zona pellucida sperm-binding protein gene (zp), egg envelope glycoprotein-like precursor, and ovarian cancer-associated gene 2 (ovca2) exhibit sex-biased expression pattern. In strategy iii, among the steroidogenic enzyme genes, steroidogenic acute regulatory protein gene (star), 17-beta-hydroxysteroid dehydrogenase type 1 gene (hsd17b1), and estradiol 17-betadehydrogenase 12 gene (hsd17b12) present sex differential expression profile. In strategy iv, cathepsins (ctss), ropporin-1-like protein gene (ropn1l), ZPA domain containing protein precursor gene (zpa), zpc5, zygote arrest protein 1 gene (zar1), wee1-like protein kinase 2 gene (wee2), P43 5S RNA-binding protein gene (42sp43), histone H2Ax gene (h2ax), and so forth display differential expression profile between male and female. The sex-related genes identified were listed in Table 2.

Sexual Dimorphic Biological Pathway.
To identify biological pathway that shows sexual dimorphism in flounder, the number of sex-biased genes was counted in different category of pathways. The result exhibits in Figure S5. It shows that the number of upregulated genes in most of metabolic pathways is much more in female than that in male. There are the most sex-biased genes in the signal transduction pathway compared with other pathways. And the significantly divergent pathways of male and female gene numbers mainly include lipid metabolism, signal transduction, translation, and cell growth death. Among them, several pathways associated with gonadal development and sex maintenance were found, such as ovarian steroidogenesis (shown as Figure 2 and Table 3), estrogen signaling pathway, progesterone-mediated oocyte maturation, prolactin signaling pathway, GnRH signaling pathway, oocyte meiosis, TGF-beta signaling pathway, steroid hormone biosynthesis, and Wnt signaling pathway.

Discussion
As a marine economic flatfish, flounder exhibit GSD + EE sex determination form, and GSD function is necessary for us to probe into. Now, little evidence can demonstrate the molecular mechanism of flounder sex determination; therefore more sex-related genes and explicit biological pathways are needed to imply the mechanism. In this study, we used RNA-seq technology to identify large quantities of sex-biased genes and illustrated its function from the perspective of biological pathways. Among these sex-biased genes, some are associated with the sex determination and gonadal development as reported. The other sex-biased genes need further study to investigate their connection with sex determination and differentiation. The sex-related genes identified in this study could provide an important clue for sex determination mechanism of flounder. In addition, the novel contigs may be from unknown gene sequence or alternatively spliced transcripts. Yano et al. [8] characterized a novel gene expressed only in testis through analyzing gonadal transcriptome of rainbow trout (O. mykiss), and further study revealed that this gene is Y chromosome sequence tightly linked with sex locus and necessary to trigger testicular differentiation. Therefore, our study also provides considerable novel sex-biased genes for further study on sex determination of flounder.

Sex Differences in Gene Expression Profiles of Gonads.
We analyzed the overall gene expression profiles of gonads and identified numerous sex-related genes. Some of the sexbiased genes are known to show sexual dimorphism between ovary and testis testifying the reliability of the selection   criterion in the results. Our analyses also found considerable previously uncharacterized sex-biased genes. Further functional characterization of these genes using transgenic overexpression, knockout strategies, and knockdown strategies may help elucidate the molecular mechanisms controlling sex determination and gonadal development in teleost.

Male-Biased Genes.
In the present study, more contigs were found to be male-biased genes than female-biased genes. For lack of genomic sequence of flounder, a large amount of male-biased contigs cannot be further assembled and annotated. Among the annotated genes, several welldocumented and important male-enriched genes are listed in more detail. Dmrt1 is expressed in the embryonic gonads of many vertebrates. It was thought to repress the female pathway through inhibition of cyp19a1a expression in Nile tilapia (O. niloticus) [34]. In medaka (O. latipes), dmy, a duplicated copy of the dmrt1 on the Y chromosome, has been confirmed as sex determining gene [5,6]. Our result shows the maledominant expression of dmrt1 in flounder as reported [25]. These reveal that dmrt1 may be an important factor for flounder testicular differentiation. Three members of sox gene family including sox9, sox8a, and sox6b were identified to be male predominantly expressed genes. In teleost, sox9 has been associated with testicular development, with sox9expressing cells differentiating into Sertoli cells [35]. Thestudy on sox8a was only reported in orange-spotted grouper (Epinephelus coioides), and it was expressed in diverse tissues and increased during testicular developmental stages [36]. No document covered the function of sox6b. AMY-1-associating protein was found to be associated with protein kinase A anchor protein 84/149 in the mitochondria in human sperm, suggesting that it plays a role in spermatogenesis [37]. Heat shock protein 90 (HSP90) interacts with steroid hormone receptors, signaling kinases, and various transcription factors [38]. However, the mechanism by which HSP90 interacts with different proteins in various pathways remains unclear.
Here in our study, we find hsp 90 genehighly expressed in both ovary and testis tissues, while hsp 90 gene exhibited upregulated expression pattern in testis. Several proteins associated with sperm mobility show absolute predominance in testis. In human, Ropporin-1-like protein is spermatogenic cell-specific protein that serves as an anchoring protein for the A-kinase anchoring protein, which may function as a regulator of both motility-and head-associated functions such as capacitation and the acrosome reaction [39]. Sperm flagellar protein 2 plays an important role in spermatogenesis and flagellar assembly. And a loss of function mutation in spef2 of mice caused the big giant head phenotype [40]. Axonemal dynein light intermediate polypeptide may take part in the formation of sperm flagella and play a dynamic role in flagellar motility. As is well known, both cytochrome c oxidase subunit I and ATP synthase F0 subunit 6 are components of the respiratory chain supplying energy for sperm mobility.

Female-Biased Genes.
Ovary is the female reproductive system and ovum-producing reproductive organ. It is responsible for the synthesis of estrogen and oogenesis.
More and more female-biased genes were found to play important roles in the above processes. In this study, nearly 66.5% of female-biased genes were annotated, of which welldocumented ovary markers such as cyp19a, foxl2, zp, cathepsins, and 42sp43 were identified as female-biased genes. Cyp19a is responsible for encoding P450 aromatase, critical enzyme catalyzing the process of transforming androgen into estrogen [41]. In teleost, cyp19a has been proved to play an important role in sex differentiation and ovarian development, and it is regarded as a reliable early marker of ovarian differentiation [42]. In addition, foxl2, encoding the activating transcription factor of cyp19a, has been identified to be uniquely expressed in female, which may be one reason for the higher expression level of cyp19a in female than in male justified by Yamaguchi et al. [24]. In mammal, zar1 gene is oocyte-specific maternal-effect gene that functions at the oocyte-to-embryo transition [43]. Wee1-like protein kinase 2 is oocyte-specific protein tyrosine kinase that phosphorylates and inhibits cyclin-dependent kinase 1 (CDK1) and acts as a key regulator of meiosis in Xenopus [44]. Histone H2Ax is reported to play an important role in chromatin remodeling and associated silencing in male mouse meiosis [45]. However, we found h2ax highly expressed in ovary rather than testes in our result, which suggests that it may act on oogenesis of female flounder. The developing oocyte is surrounded by an acellular envelope that is composed of zona pellucida proteins. It is reported that zpa, zpc, and zp3 genes and these proteins participate in taxon-specific spermegg binding during fertilization process and protect embryo at early developmental stage in mammals [46]. Zps were also identified as female-biased genes in this study. Vitellogenesis is the principal event responsible for the enormous growth of oocytes in many teleosts, during which most nutritive products are taken up and stored for developing embryo. Vitellogenin receptor is involved in uptake of vitellogenin by endocytosis. Enzymes such as cathepsins are responsible for the degradation of vitellogenin into yolk protein for storage in the oocyte [47]. P43 5S RNA-binding protein is combined with 5S rRNA to comprise 42S ribonucleoprotein storage particle. In addition, transcription factor IIIA acts as both a positive transcription factor for 5S RNA genes and a specific RNA-binding protein that complex with 5S RNA in oocytes to form the 7S ribonucleoprotein storage particle [48]. According to Diaz De Cerio et al., 5 S RNA and associate protein could constitute a sensitive and universal marker of oogenesis and oocyte differentiation in fish [49].

Sex Steroids and Biosynthetic Pathway.
Gonad is an important organ that is responsible for producing gametes and sex hormones. Steroid hormones are small, hydrophobic hormones that can permeate membranes and therefore can bind to specific nuclear receptors, including estrogen receptors (ER) and androgen [50]. It can form a complex with a hormone receptor and enter the nucleus where this complex can bind to DNA and result in the translation of specific mRNA and proteins. These can give rise to specific physiological responses including sex differentiation and germ-cell development [51]. In the present study, 20 genes were identified to be involved in flounder ovarian steroidogenesis pathway. In this pathway, all steroids are synthesized from cholesterol, whose transportation from intracellular source into the mitochondria is the rate-limiting step in steroidogenesis [52]. This transmembrane transport process is facilitated by the steroidogenic acute regulatory protein [53]. And the present study identified star as malebiased gene. Male hormone testosterone is formed from pregnenolone by two pathways, delta5 pathway via dehydroepiandrosterone and delta4 pathway via androstenedione. The enzyme P450c17 is responsible for the 17,20-lyase and 17-alphahydroxylase activities in respective pathways. Two forms of P450c17 were identified in some teleost, which were P450c17-I and P450c17-II, encoded by cyp17-I and cyp17-II genes, respectively [54]. P450c17-I possesses both hydroxylase and lyase activities, while P450c17-II only has hydroxylase activity. In medaka, P450c17-I is essential for the production of estradiol-17 (E2) during oocyte growth, while P450c17-II plays a vital role in production of 17 , 20 -dihydroxy-4-pregnen-3-one (17 , 20 -DP) during oocyte maturation [55]. In our study, the transcripts of cyp17-II were only found. According to Ding et al., the variation trends of T and E2 level were consistent with the cyp17-II expression pattern in flounder ovary [26]. Nevertheless, cyp19a still plays a central role in the synthesis of E2. The suppression of cyp19a expression at male-inducing temperature leads to masculinization, which further supports the importance of this enzyme and its product (estrogen) in flounder ovarian differentiation [23]. Additionally, 17-beta-hydroxysteroid dehydrogenase (hsd17b) and 3-beta-hydroxysteroid dehydrogenase (hsd3b) were identified to be female predominantly expressed genes, which indicate their important role in synthesis of estrogen.
The steroidogenic enzymes expression profile may help us to demonstrate the difference of steroid level between male and female. Based on our previous study, both testosterone and estradiol-17 (E2) levels are on rise from ovarian developmental stage I to stage IV. The decrease of E2 level is detected in temperature-induced masculinized groups compared to control groups [56]. During testicular developmental stage, E2 in the serum stay in low level while T level varies significantly during different developmental stage. The variation of sex steroids level may be tightly linked with gonadal development and maturation of germ cells.

Conclusion
This is the first report of flounder gonadal transcriptome using RNA-seq technology. We generated a large number of ESTs collection and identified numerous differentially expressed genes between ovary and testis. According to annotation information, sex-related biological pathways including ovarian steroidogenesis were found. The dimorphic expression of overall sex-related genes provides further insights into sexual difference and gonadal development. Our result also provides an archive for further study on molecular mechanism underlying sex determination.