BioMed Research International

Volume 2014 (2014), Article ID 291067, 10 pages

http://dx.doi.org/10.1155/2014/291067

Research Article

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

Zhaofei Fan,1,2 Feng You,¹ Lijuan Wang,1,2 Shenda Weng,¹ Zhihao Wu,¹ Jinwei Hu,1,2 Yuxia Zou,¹ Xungang Tan,¹ and Peijun Zhang¹

¹Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences, 7 Nanhai Road, Qingdao 266071, China

²University of the Chinese Academy of Sciences, Beijing 10049, China

Received 7 April 2014; Revised 12 June 2014; Accepted 15 June 2014; Published 6 July 2014

Academic Editor: Tatsuya Akutsu

Copyright © 2014 Zhaofei Fan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. 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]. Sex-related 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 sex-related 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.

2. Materials and Methods

2.1. Fish

Adult and juvenile flounder (1240 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³ 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 (4050 cm in TL) and 3 mature female flounder (5060 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.

2.2. 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].

2.3. 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 first-strand 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₄ (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.

2.4. 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 -value of top hit lower than 1 . Furthermore, GoPipe software was used to perform blastp (cut-off -value of 1 ) 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 -value of 1 ) 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.

2.5. 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, “’’ and “|log2 (RPKM_XX/RPKM_XY)| ’’ 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 -value of 1 ) 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 sex-biased genes.

3. Results

3.1. Assembly, Annotation, and Bioinformatical Analysis

Approximately 20.1 million and 22.4 million reads were obtained from male library and female library, respectively. After quality filtering, about 98.1% read of male and 99.0% read of female remain to be qualified for assembling. Sequencing saturation distribution (see Figure S1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/291067) and genes coverage statistic (Figure S2) analysis justified a deep sequencing coverage sufficient for the quantitative analysis of gene expression profiles. The clean reads were jointly assembled into 97,233 contigs with N50 of 809 bp and average length of 603 bp. The size distribution of assembled contigs was presented in Figure 1. Mapped to protein database, nearly 22.31% (21,697) of contigs were matched with known existing protein (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. About 21 categories of biological process were assigned for 49,204 contigs, 11 categories of cellular component were assigned for 49,253 contigs, and 25 categories of molecular function were assigned for 41,788 contigs. Most of the cellular component genes were associated with cells and intracellular components, and most of the molecular function genes were associated with binding and catalytic activity. Mapped to reference canonical pathways in the KEGG database, 11,936 contigs were assigned with KO numbers. Through the KO number of the predicted protein, 328 metabolic pathways were constructed with various degrees (Table S1_KEGG pathway). Figure S4 showed the category distribution of biological pathways, in which the most enriched pathway is associated with signal transduction.

Figure 1: Length distribution of contigs.

3.2. 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-beta-dehydrogenase 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.

Table 1: Overview of transcriptome analysis.

Total number of reads
42,640,333

    Number of male reads
20,161,959

    Number of female reads
22,478,374

    Number of male reads in contigs
14,266,081

    Number of female reads in contigs
20,939,047

Number of contigs from assembly
97,233

    Number of annonated contigs
21,697

    Number of contigs >400 bp
41,135

    Number of male only contigs
25,226

    Number of female only contigs
10,750

    Number of mixed contigs
61,257

Sex-biased genes with -value ≤ 1 × 10⁻³ AND ratio ≥ 2
13,644

    Male-biased contigs
10,348

    Annonated male-biased genes
887

    Female-biased contigs
3296

    Annonated female-biased genes
2193

Table 2: Representative sex-biased genes involved in sex determination and development.

Contigs
Length
Function
(fold)
-value
Sex
Gene

Comp134314_c0_seq1
2117
P450 aromatase
8.1274

Female
cyp19a

Comp128200_c0_seq1
1721
Forkhead transcription factor L2
7.7734

Female
foxl2

Comp129336_c0_seq1
697
Orphan nuclear receptor Dax1
2.3341

Female
dax1

Comp135979_c0_seq1
3970
Vitellogenin receptor
2.4058
0
Female
vtgr

Comp131376_c0_seq1
1862
Zona pellucida sperm-binding protein 4
6.733436
0
Female
zp4

Comp129690_c1_seq1
1123
Zona pellucid sperm-binding protein 3
6.8619
0
Female
zp3

Comp126058_c0_seq1
1449
Zona pellucida sperm-binding protein
6.1157
0
Female
zp

Comp96429_c0_seq1
668
Histone H2A.x
6.3287
0
Female
h2a.x

Comp124668_c0_seq1
1372
Zygote arrest protein 1
6.556118
0
Female
zar1

Comp138502_c2_seq1
2515
ZPA domain containing protein precursor
7.7351
0
Female
zpa

Comp131365_c0_seq1
2442
ZPC5
6.5016
0
Female
zpc5

Comp130117_c0_seq1
1092
High choriolytic enzyme 2
7.5466
0
Female
hce2

Comp135460_c1_seq1
4433
Protein fem-1 homolog C
2.2905

Female
fem1c

Comp130482_c0_seq1
1936
Frizzled-3
3.2588

Female
fzd3

Comp126332_c0_seq1
1789
P43 5S RNA-binding protein
8.3272
0
Female
42sp43

Comp134412_c0_seq1
3026
Transcription factor IIIA
3.7385
0
Female
gtf3a

Comp129656_c0_seq1
2020
Wee1-like protein kinase 2
6.9006
0
Female
wee2

Comp126733_c0_seq1
1431
Cathepsin S precursor
5.1171
0
Female
ctss

Comp137073_c0_seq1
2247
G2/mitotic-specific cyclin-B1-like
6.7561
0
Female
ccnb1

Comp126534_c0_seq1
1128
Ovarian cancer-associated gene 2 protein
2.744

Female
ovca2

Comp130477_c0_seq1
1312
17-beta-Hydroxysteroid

dehydrogenase type 1
5.2759

Female
hsd17b1

Comp128629_c0_seq1
1286
SRY-box containing protein 9
−4.0548

Male
sox9

Comp138102_c0_seq1
1849
Mullerian inhibiting substance
−5.7273

Male
mis

Comp135961_c0_seq1
2582
Sox8a
−4.3658
0
Male
sox8a

Comp125855_c0_seq1
1045
Steroidogenic acute regulatory protein
−2.1114

Male
star

Comp70046_c1_seq1
523
SRY-box containing gene 6b
−2.9959

Male
sox6b

Comp125855_c0_seq1
1045
Steroidogenic acute regulatory protein
−2.1114

Male
star

Comp70129_c0_seq1
1376
Cytochrome c oxidase subunit I
−13.8122
0
Male
cox I

Comp127682_c0_seq1
960
ATP synthase F0 subunit 6
−12.5508

Male
atp5g

Comp127263_c1_seq1
1000
Heat shock protein 90 alpha
−6.5865

Male
hsp90α

Comp132152_c1_seq1
3496
Sperm flagellar protein 2
−3.2844

Male
spef

Comp138413_c1_seq1
2470
AMY-1-associating protein
−2.2291

Male
aat1

Comp131350_c0_seq1
1398
Axonemal dynein light intermediate polypeptide 1
−9.422
0
Male
dnali1

Comp128419_c1_seq1
659
Ropporin-1-like protein
−8.2764

Male
ropn1l

Note: Fold indicates RPKM (female)/RPKM (male); -value indicates false discove rate (FDR).

3.3. 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.

Table 3: Genes identified to take part in the ovarian steroidogenesis pathway.

Contigs
Gene function
Protein
-value
(fold)

Comp1153048_c0_seq1
Luteinizing hormone receptor
LHR

NA/female

Comp67657_c0_seq1
Guanine nucleotide-binding protein G(s) subunit alpha
Gs

0.299

Comp127442_c3_seq1
Insulin receptor
INSR

2.911

Comp123852_c0_seq1
Low-density lipoprotein receptor-related protein 1
LDLR
0
1.607

Comp137524_c1_seq1
Follicle stimulating hormone receptor II
FSHR
0
−1.703

Comp108327_c0_seq1
cAMP-dependent protein kinase catalytic subunit PRKX
PKA

1.911

Comp138426_c0_seq1
Acyl-protein thioesterase 2
ACOT2

1.514

Comp1185767_c0_seq1
Cholesterol side chain cleavage cytochrome P450
CYP11A1

NA/male

Comp945537_c0_seq1
17-alpha-Hydroxylase
CYP17

1.446

Comp125855_c0_seq1
Steroidogenic acute regulatory protein, mitochondrial
StAR

−2.111

Comp134314_c0_seq1
P450 aromatase
CYP19A1
0
8.127

Comp123514_c0_seq1
3-beta-Hydroxysteroid dehydrogenase type II
HSD3b

2.288

Comp134625_c0_seq1
Estradiol 17-beta-dehydrogenase 12-A
HSD17b

4.625

Comp131665_c2_seq1
Cytochrome P450 1B1
CYP1B1
0
−1.964

Comp123774_c0_seq1
Cytochrome P450 1A
CYP1A1
0
−2.691

Comp126639_c0_seq1
Catechol O-methyltransferase
COMT

−1.249

Comp507499_c0_seq1
Cytosolic phospholipase A2
CPLA2

−2.548

Comp80365_c0_seq1
Cytochrome P450 2J2
CYP2J

0.451

Comp133304_c0_seq1
Scavenger receptor class B member
SCARB1
0
0.919

Comp136599_c2_seq1
Adenylate cyclase type 1
ADCY1

4.797

Note: NA indicates solely identification in female or male; fold indicates RPKM (female)/RPKM (male).

Figure 2: Ovarian steroidogenesis pathway. Green background indicates genes identified in flounder; red background indicates female-biased genes; blue background indicates male-biased genes.

4. 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.

4.1. 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 sex-biased 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.

4.1.1. 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 well-documented 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 male-dominant 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 sox9-expressing 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.

4.1.2. 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 well-documented 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 sperm-egg 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].

4.2. 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 male-biased 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-alpha-hydroxylase 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.

5. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This study was supported by the National High Technology Research and Development Program of China (863 Programme, no. 2012AA092203), the National Natural Science Foundation of China (no. 41276171), the National Key Basic Program of Science and Technology-Platforms of Aquaculture Stock Resources, and the National Flatfish Industry System Construction Programme (no. nycytx-50-G03).

References

G. S. Helfman, B. B. Collette, D. E. Facey, and B. Bowen, The Diversity of Fishes. Biology, Evolution, and Ecology, Blackwell Science, Malden, Mass, USA, 2nd edition, 2009.

J. A. Luckenbach, R. J. Borski, H. V. Daniels, and J. Godwin, “Sex determination in flatfishes: mechanisms and environmental influences,” Seminars in Cell & Developmental Biology, vol. 20, no. 3, pp. 256–263, 2009.

A. Cnaani and B. Levavi-Sivan, “Sexual development in fish, practical applications for aquaculture,” Sexual Development, vol. 3, no. 2-3, pp. 164–175, 2009.

N. Ospina-Álvarez and F. Piferrer, “Temperature-dependent sex determination in fish revisited: Prevalence, a single sex ratio response pattern, and possible effects of climate change,” PLoS ONE, vol. 3, no. 7, Article ID e2837, 2008.

I. Nanda, M. Kondo, U. Hornung et al., “A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes,” Proceedings of the National Academy of Sciences of the United States of America, vol. 99, no. 18, pp. 11778–11783, 2002.

M. Matsuda, Y. Nagahama, A. Shinomiya et al., “DMY is a Y-specific DM-domain gene required for male development in the medaka fish,” Nature, vol. 417, no. 6888, pp. 559–563, 2002.

R. S. Hattori, Y. Murai, M. Oura et al., “A Y-linked anti-Müllerian hormone duplication takes over a critical role in sex determination,” Proceedings of the National Academy of Sciences of the United States of America, vol. 109, no. 8, pp. 2955–2959, 2012.

A. Yano, R. Guyomard, B. Nicol et al., “An immune-related gene evolved into the master sex-determining gene in rainbow trout, Oncorhynchus mykiss,” Current Biology, vol. 22, no. 15, pp. 1423–1428, 2012.

T. Myosho, H. Otake, H. Masuyama et al., “Tracing the emergence of a novel sex-determining gene in medaka, Oryzias luzonensis,” Genetics, vol. 191, no. 1, pp. 163–170, 2012.

T. Kamiya, W. Kai, S. Tasumi et al., “A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger Pufferfish, Takifugu rubripes (Fugu),” PLoS Genetics, vol. 8, no. 7, Article ID e1002798, 2012.

S. Chen, G. Zhang, C. Shao, Q. Huang, and G. Liu, “Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle,” Nature Genetics, vol. 46, pp. 253–260, 2014.

G. Sandra and M. Norma, “Sexual determination and differentiation in teleost fish,” Reviews in Fish Biology and Fisheries, vol. 20, no. 1, pp. 101–121, 2010.

A. Cutting, J. Chue, and C. A. Smith, “Just how conserved is vertebrate sex determination?” Developmental Dynamics, vol. 242, no. 4, pp. 380–387, 2013.

Z. Zhang, Y. Wang, S. Wang et al., “Transcriptome analysis of female and male Xiphophorus maculatus Jp 163 A,” PLoS ONE, vol. 6, no. 4, Article ID e18379, 2011.

M. Salem, C. E. Rexroad III, J. Wang, G. H. Thorgaard, and J. Yao, “Characterization of the rainbow trout transcriptome using Sanger and 454-pyrosequencing approaches,” BMC Genomics, vol. 11, no. 1, article 564, 2010.

W. Tao, J. Yuan, L. Zhou et al., “Characterization of gonadal transcriptomes from Nile Tilapia (Oreochromis niloticus) reveals differentially expressed genes,” PLoS ONE, vol. 8, no. 5, Article ID e63604, 2013.

L. Sun, C. Wang, L. Huang, M. Wu, and Z. Zuo, “Transcriptome analysis of male and female Sebastiscus marmoratus,” PLoS ONE, vol. 7, no. 11, Article ID e50676, 2012.

F. Sun, S. Liu, X. Gao et al., “Male-biased genes in catfish as revealed by RNA- seq analysis of the testis transriptome,” PLoS ONE, vol. 8, no. 7, Article ID e68452, 2013.

L. Ribas, B. G. Pardo, C. Fernández et al., “A combined strategy involving Sanger and 454 pyrosequencing increases genomic resources to aid in the management of reproduction, disease control and genetic selection in the turbot (Scophthalmus maximus),” BMC Genomics, vol. 14, no. 1, article 180, 2013.

E. Yamamoto, “Studies on sex-manipulation and production of cloned populations in hirame, Paralichthys olivaceus,” Bulletin of the Tottori Prefectual Fisheries Experimental Station, vol. 34, pp. 1–145, 1995.

E. Yamamoto, “Studies on sex-manipulation and production of cloned populations in hirame, Paralichthys olivaceus (Temminck et Schlegel),” Aquaculture, vol. 173, no. 1–4, pp. 235–246, 1999.

P. Sun, F. You, M. Liu et al., “Steroid sex hormone dynamics during estradiol-17β induced gonadal differentiation in Paralichthys olivaceus (Teleostei),” Chinese Journal of Oceanology and Limnology, vol. 28, no. 2, pp. 254–259, 2010.

T. Kitano, K. Takamune, T. Kobayashi, Y. Nagahama, and S.-I. Abe, “Suppression of P450 aromatase gene expression in sex-reversed males produced by rearing genetically female larvae at a high water temperature during a period of sex differentiation in the Japanese flounder (Paralichthys olivaceus),” Journal of Molecular Endocrinology, vol. 23, no. 2, pp. 167–176, 1999.

T. Yamaguchi, S. Yamaguchi, T. Hirai, and T. Kitano, “Follicle-stimulating hormone signaling and Foxl2 are involved in transcriptional regulation of aromatase gene during gonadal sex differentiation in Japanese flounder, Paralichthys olivaceus,” Biochemical and Biophysical Research Communications, vol. 359, no. 4, pp. 935–940, 2007.

A. Y. Wen, F. You, and P. Sun, “CpG methylation of dmrt1 and cyp19a promoters in relation to their sexual dimorphic expression in the Japanese flounder Paralichthys olivaceus,” Journal of Fish Biology, vol. 84, no. 1, pp. 193–205, 2014.

Y. Ding, F. He, H. Wen et al., “DNA methylation status of cyp17-II gene correlated with its expression pattern and reproductive endocrinology during ovarian development stages of Japanese flounder (Paralichthys olivaceus),” Gene, vol. 527, no. 1, pp. 82–88, 2013.

A. Wen, F. You, X. Tan et al., “Expression pattern of dmrt4 from olive flounder (Paralichthys olivaceus) in adult gonads and during embryogenesis,” Fish Physiology and Biochemistry, vol. 35, no. 3, pp. 421–433, 2009.

P. Sun, F. You, L. Zhang, et al., “Histological evaluation of gonadal differentiation in olive flouder (Paralichthys olivaceus),” Marine Sciences, vol. 33, pp. 53–58, 2009.

M. Radonic and G. J. Macchi, “Gonadal sex differentiation in cultured juvenile flounder, Paralichthys orbignyanus (Valenciennes, 1839),” Journal of the World Aquaculture Society, vol. 40, no. 1, pp. 129–133, 2009.

Z. Z. Chen, C. H. Xue, S. Zhu et al., “GoPipe: streamlined gene ontology annotation for batch anonymous sequences with statistics,” Progress in Biochemistry and Biophysics, vol. 32, no. 2, pp. 187–191, 2005.

M. Kanehisa, S. Goto, M. Furumichi, M. Tanabe, and M. Hirakawa, “KEGG for representation and analysis of molecular networks involving diseases and drugs,” Nucleic Acids Research, vol. 38, no. 1, pp. D355–D360, 2009.

A. Mortazavi, B. A. Williams, K. McCue, L. Schaeffer, and B. Wold, “Mapping and quantifying mammalian transcriptomes by RNA-Seq,” Nature Methods, vol. 5, no. 7, pp. 621–628, 2008.

L. Wang, Z. Feng, X. Wang, and X. Zhang, “DEGseq: an R package for identifying differentially expressed genes from RNA-seq data,” Bioinformatics, vol. 26, no. 1, pp. 136–138, 2010.

D. S. Wang, L. Y. Zhou, T. Kobayashi et al., “Doublesex- and Mab-3-related transcription factor-1 repression of aromatase transcription, a possible mechanism favoring the male pathway in tilapia,” Endocrinology, vol. 151, no. 3, pp. 1331–1340, 2010.

J. Kent, S. C. Wheatley, J. E. Andrews, A. H. Sinclair, and P. Koopman, “A male-specific role for SOX9 in vertebrate sex determination,” Development, vol. 122, no. 9, pp. 2813–2822, 1996.

Q. Liu, H. Lu, L. Zhang, J. Xie, W. Shen, and W. Zhang, “Homologues of sox8 and sox10 in the orange-spotted grouper Epinephelus coioides: sequences, expression patterns, and their effects on cyp19a1a promoter activities in vitro,” Comparative Biochemistry and Physiology B, vol. 163, no. 1, pp. 86–95, 2012.

H. Yukitake, M. Furusawa, T. Taira, S. M. M. Iguchi-Ariga, and H. Ariga, “AAT-1, a novel testis-specific AMY-1-binding protein, forms a quaternary complex with AMY-1, a-kinase anchor protein 84, and a regulatory subunit of cAMP-dependent protein kinase and is phosphorylated by its kinase,” Journal of Biological Chemistry, vol. 277, no. 47, pp. 45480–45492, 2002.

W. Liu, F. X. Zhang, M. J. Cai et al., “The hormone-dependent function of Hsp90 in the crosstalk between 20-hydroxyecdysone and juvenile hormone signaling pathways in insects is determined by differential phosphorylation and protein interactions,” Biochimica et Biophysica Acta, vol. 1830, no. 11, pp. 5184–5192, 2013.

D. W. Carr, A. Fujita, C. L. Stentz, G. A. Liberty, G. E. Olson, and S. Narumiya, “Identification of sperm-specific proteins that interact with A-kinase anchoring proteins in a manner similar to the type II regulatory subunit of PKA,” Journal of Biological Chemistry, vol. 276, no. 20, pp. 17332–17338, 2001.

A. Sironen, J. Hansen, B. Thomsen et al., “Expression of SPEF2 during mouse spermatogenesis and identification of IFT20 as an interacting protein,” Biology of Reproduction, vol. 82, no. 3, pp. 580–590, 2010.

A. Conley and M. Hinshelwood, “Mammalian aromatases,” Reproduction, vol. 121, no. 5, pp. 685–695, 2001.

Y. Guiguen, A. Fostier, F. Piferrer, and C.-F. Chang, “Ovarian aromatase and estrogens: a pivotal role for gonadal sex differentiation and sex change in fish,” General and Comparative Endocrinology, vol. 165, no. 3, pp. 352–366, 2010.

X. Wu, M. M. Viveiros, J. J. Eppig, Y. Bai, S. L. Fitzpatrick, and M. M. Matzuk, “Zygote arrest 1 (Zar1) is a novel maternal-effect gene critical for the oocyte-to-embryo transition,” Nature Genetics, vol. 33, no. 2, pp. 187–191, 2003.

K. Okamoto, N. Nakajo, and N. Sagata, “The existence of two distinct Wee1 isoforms in Xenopus: implications for the developmental regulation of the cell cycle,” EMBO Journal, vol. 21, no. 10, pp. 2472–2484, 2002.

O. Fernandez-Capetillo, S. K. Mahadevaiah, A. Celeste et al., “H2AX is required for chromatin remodeling and inactivation of sex chromosomes in male mouse meiosis,” Developmental Cell, vol. 4, no. 4, pp. 497–508, 2003.

C. Modig, T. Modesto, A. Canario, J. Cerdà, J. von Hofsten, and P. Olsson, “Molecular characterization and expression pattern of zona pellucida proteins in gilthead seabream (Sparus aurata),” Biology of Reproduction, vol. 75, no. 5, pp. 717–725, 2006.

S. Brooks, C. R. Tyler, and J. P. Sumpter, “Egg quality in fish: what makes a good egg?” Reviews in Fish Biology and Fisheries, vol. 7, no. 4, pp. 387–416, 1997.

H. Schneider, M. Dabauvalle, N. Wilken, and U. Scheer, “Visualizing protein interactions involved in the formation of the 42S RNP storage particle of Xenopus oocytes,” Biology of the Cell, vol. 102, no. 8, pp. 469–478, 2010.

O. Diaz De Cerio, I. Rojo-Bartolomé, C. Bizarro, M. Ortiz-Zarragoitia, and I. Cancio, “5S rRNA and accompanying proteins in gonads: powerful markers to identify sex and reproductive endocrine disruption in fish,” Environmental Science and Technology, vol. 46, no. 14, pp. 7763–7771, 2012.

G. E. Sandra and M. M. Norma, “Sexual determination and differentiation in teleost fish,” Reviews in Fish Biology and Fisheries, vol. 20, no. 1, pp. 101–121, 2010.

R. H. Devlin and Y. Nagahama, “Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences,” Aquaculture, vol. 208, no. 3-4, pp. 191–364, 2002.

M. B. Rone, J. Fan, and V. Papadopoulos, “Cholesterol transport in steroid biosynthesis: Role of protein-protein interactions and implications in disease states,” Biochimica et Biophysica Acta, vol. 1791, no. 7, pp. 646–658, 2009.

P. R. Manna, M. T. Dyson, and D. M. Stocco, “Regulation of the steroidogenic acute regulatory protein gene expression: present and future perspectives,” Molecular Human Reproduction, vol. 15, no. 6, pp. 321–333, 2009.

L.-Y. Zhou, D.-S. Wang, T. Kobayashi et al., “A novel type of P450c17 lacking the lyase activity is responsible for C21-steroid biosynthesis in the fish ovary and head kidney,” Endocrinology, vol. 148, no. 9, pp. 4282–4291, 2007.

L. Y. Zhou, D. S. Wang, Y. Shibata, B. Paul-Prasanth, A. Suzuki, and Y. Nagahama, “Characterization, expression and transcriptional regulation of P450c17-I and -II in the medaka, Oryzias latipes,” Biochemical and Biophysical Research Communications, vol. 362, no. 3, pp. 619–625, 2007.

P. Sun, F. You, D. Ma, J. Li, and P. Zhang, “Sex steroid changes during temperature-induced gonadal differentiation in Paralichthys olivaceus (Temminck & Schegel, 1846),” Journal of Applied Ichthyology, vol. 29, no. 4, pp. 886–890, 2013.