Next Article in Journal
Multifaceted Roles of ALG-2 in Ca2+-Regulated Membrane Trafficking
Next Article in Special Issue
Transcriptome Sequencing and De Novo Assembly of Golden Cuttlefish Sepia esculenta Hoyle
Previous Article in Journal
Aquaporin-4 in Astroglial Cells in the CNS and Supporting Cells of Sensory Organs—A Comparative Perspective
Previous Article in Special Issue
Age-Specific Lipid and Fatty Acid Profiles of Atlantic Salmon Juveniles in the Varzuga River
Article

Sexually Dimorphic Gene Expression Associated with Growth and Reproduction of Tongue Sole (Cynoglossus semilaevis) Revealed by Brain Transcriptome Analysis

1
Heilongjiang River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Harbin 150070, China
2
Environmental Engineering Program, Department of Civil Engineering, Auburn University, Auburn, AL 36849, USA
3
School of Computer Science and Technology, Harbin Institute of Technology, Harbin 150001, China
*
Author to whom correspondence should be addressed.
Academic Editors: Jun Li and Li Lin
Int. J. Mol. Sci. 2016, 17(9), 1402; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms17091402
Received: 25 March 2016 / Revised: 5 August 2016 / Accepted: 19 August 2016 / Published: 26 August 2016
(This article belongs to the Special Issue Fish Molecular Biology)

Abstract

In this study, we performed a comprehensive analysis of the transcriptome of one- and two-year-old male and female brains of Cynoglossus semilaevis by high-throughput Illumina sequencing. A total of 77,066 transcripts, corresponding to 21,475 unigenes, were obtained with a N50 value of 4349 bp. Of these unigenes, 33 genes were found to have significant differential expression and potentially associated with growth, from which 18 genes were down-regulated and 12 genes were up-regulated in two-year-old males, most of these genes had no significant differences in expression among one-year-old males and females and two-year-old females. A similar analysis was conducted to look for genes associated with reproduction; 25 genes were identified, among them, five genes were found to be down regulated and 20 genes up regulated in two-year-old males, again, most of the genes had no significant expression differences among the other three. The performance of up regulated genes in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was significantly different between two-year-old males and females. Males had a high gene expression in genetic information processing, while female’s highly expressed genes were mainly enriched on organismal systems. Our work identified a set of sex-biased genes potentially associated with growth and reproduction that might be the candidate factors affecting sexual dimorphism of tongue sole, laying the foundation to understand the complex process of sex determination of this economic valuable species.
Keywords: sexually dimorphic; tongue sole; transcriptome analysis sexually dimorphic; tongue sole; transcriptome analysis

1. Introduction

Tongue sole (Cynoglossus semilaevis) is a valuable marine fish in China, which is receiving more attention because of its good taste and nutritional value. However, the species is decreasing in number in the wild despite increasing market demand, which favors more production of the species [1]. Apart from its increasing economic value, it also has an intriguing reproductive biology. Genetically, tongue sole has ZW sex-determination system [2]. Males are the homogametic sex (ZZ), while females are the heterogametic sex (ZW). Male and female tongue sole are considerably different in size and growth rate: mature females are twice larger in length and six times greater in weight than their male counterparts [3]. Thus, understanding the underpinning of sexual dimorphism and sex determination in this species is essential for developing methods to boost its productivity to meet the aquaculture market demands.
Although sexual determination in some fish is partly determined by the environment [4,5,6], brain is still a major organ involved in fish growth and reproduction [7], and plays a key role in sexual dimorphism and sex determination, including the regulation of reproduction, maturation and sexual behavior in both sexes. The influence of brain on development is mainly mediated by the Brain–Pituitary–Gonadal (BPG) axis [8]. Gahr [9] found that male quail transplanted with female forebrain do not develop normal testes, indicating that a male brain is necessary for normal testes development. This demonstrates that sexual dimorphism in brain is not entirely determined by hormones secreted by the gonads, and brain develops differently between sexes even before the gonad development [10,11,12]. Therefore, understanding the mechanism of sexual dimorphism in brain is fundamental for enhancing our knowledge of development processes and exploring the genetic toolkit to maintain or change the difference between sexes.
Yang and colleagues [13] found that relatively small changes in gene expression might control functional sexual dimorphism in the brain, suggesting that differences in a few of those genes may become biological significant when interact with other genes in functional gene networks. Thus, uncovering the sex-biased genes are thought to be necessary when studying the mechanism of sexual dimorphism, and these genes may evolve differently under selection pressure, in turn, acting on the sexual phenotype [14].
For researchers, the major challenge is to explain the huge list of differentially expressed genes (DEGs) under certain conditions. A decade ago, automatic functional classification method using Gene Ontology (GO) was proposed to help study the function of these DEGs [15]. The first step is to group the list of DEGs into functional groups, which can offer insight into the biological mechanisms at the given condition. In addition, a well-characterized transcriptome is considered to be helpful to identify functional genes underlying the certain phenotypic variation in several ways.
Gonads, as the main sex organ, are often utilized to get a comprehensive overview of differentially expressed genes responsible for sexual dimorphism. Up to now, thousands of sex-biased genes have been identified from gonads in many species, including fruit fly [16], mouse [17] and some teleost, such as zebrafish [18], Nile tilapia [19], sharpsnout seabream [20], guppy [21], and yellow catfish [22]. However, sex-biased genes expressed in brain have been less well studied.
For the sexually dimorphic development of tongue sole, studies have mainly focused on the hormones expressed in brain known to influence the BPG-axis such as pituitary adenylate cyclase-activating polypeptide (PACAP) [23,24], growth hormone-releasing hormone (GHRH) [23,24] and growth hormone (GH) [25,26]. However, few studies have been published at the transcriptomic level to identify and characterize sex-biased genes related to growth and reproduction in tongue sole brain. Recently, the whole genome of tongue sole was published [2], enabling whole genome-based transcriptomic analyses.
In this study, we seek to identify the molecular mechanism underlying sexual dimorphism of the brain transcriptome of tongue sole. Fish used for this project were one- and two-year-old tongue sole. The average body weights of one-year-old male and female were 54 and 99 g, and the average body lengths were 21 and 25 cm, respectively. For two-year-old fish, the average weights were 337 and 1380 g, and average lengths were 27 and 42 cm, for male and female respectively. Females were much bigger than males. Two-year-old females were still immature, while males were mature. We used the RNA-Sequencing approach [27] aimed at identifying the differentially expressed genes of tongue sole. Brain transcriptomic profiles of one- and two-year-old males and females were analyzed and revealed a number of sexually differentially expressed genes potentially associated with growth and reproduction. These data will hopefully provide a useful genetic resource for further sex determination studies on tongue sole.

2. Results

2.1. Transcriptome Sequencing and Assembly

To obtain comprehensive information of gene expression between sexes, we sequenced four transcriptomes derived from female and male brains at one and two-year-old tongue sole by using Illumina Hiseq2000 technology. A workflow with the transcriptome data analysis procedures is shown in Figure 1.
In total, 201,822,072 paired-end reads were generated, encompassing 20.2 Gb of sequence. Low quality reads and reads with ambiguous nucleotides were discarded, leaving 191,631,250 clean reads for transcriptome assembly and analysis. To assess our transcriptome assembly, we mapped our clean paired-end reads to the C. semilaevis whole genome. Approximately 84.2% of them exhibited significant hits to the genome. These reads were assembled into 77,066 transcripts in total, corresponding to 21,475 unigenes, with an average length of 3059 bp, the N50 value of 4349 bp and the N90 value of 1638 bp. Transcriptome data obtained from four samples has been uploaded to NCBI SRA site, with accession numbers of SRR2046146, SRR2046147, SRR2046153, and SRR2046154, respectively.

2.2. Gene Identification and Functional Annotation

All transcribed genes were functionally annotated by blasting the NCBI non-redundant database. A total of 10,185 genes were annotated with at least one significant hit (E value < 1 × 10−5). Of these genes, 5467 genes were further annotated by Gene Ontology (GO) with 5845 GO terms. A histogram of the number and percent of genes falling into the GO categories is given in Figure S1. The category of biological process contained the most genes, followed by molecular function and cellular component. Briefly, in biological process, cellular processes and metabolic processes were the most represented categories; with a total of 3553 genes assigned to cellular process (GO: 0009987, 64.6%) and 2756 genes to metabolic process (GO: 0008152, 50.2%). In the category with cellular component, cell (GO: 0005623, 59.2%) and organelle (GO: 0043226, 30.5%) were annotated with most abundant genes. Finally, in the category of molecular functions, binding and catalytic activity were two most frequent classes.
Functional pathway analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was also conducted as a complementary approach to categorize and annotate transcribed genes. Of the 21,475 transcribed genes, 7879 genes were assigned to KEGG orthologs, and then categorized into six functional groups (Table S1A). In short, of the genes with KEGG annotation, 1929 genes (22.5%) were assigned to metabolism, 1394 (17.7%) to genetic information processing, and 2774 genes (35.2%) were classified to environmental information processing. Cellular processes, organismal systems and human diseases contained 1883 (23.8%), 2944 (37.3%), and 2695 (34.2%) KEGG annotated genes, respectively.

2.3. Differentially Expressed Genes (DEGs) and Gene Enrichment Analysis

Gene expression profiling of the male vs. female brain at two developmental stages are plotted in Figure 2. From the total of 21,475 transcribed genes, 771 genes were found to have sex-specific expression pattern in one-year-old fish, 469 of them were up-regulated in the male and the rest in the female. For the two-year-old fish, the number of DEGs (5498 genes) was about seven times greater than that in one-year-old fish. Among those DEGs, 3025 genes were found to be up-regulated in male, and the rest in female. More up-regulated genes were found in M2 than F2, which might be related to the maturation of males in this stage. Overall, co-expressed genes (CEGs) still made up the majority of all genes at both developmental stages. Transcriptome profiling and differentially expressed genes were listed in supplementary Table S2. DEGs, CEGs and NEGs between male and female tongue sole were shown in supplementary Table S4.
Gene enrichment analysis was conducted to further understand the biological functions of the DEGs. All DEGs were mapped to the terms of Gene Ontology and compared with the background of the whole transcriptome (Table 1). More DEGs in two-year-old fish were identified (711 up-regulated genes for male (M2), and 619 genes for female (F2)) involved in Gene Ontology enrichment analysis comparing to one-year-old fish (104 genes for male (M1), and 63 genes for female (F1)). DEGs were significantly enriched in several GO terms. In biological process, two-year-old male tongue sole (M2) had a strikingly greater gene expression in reproduction (20 genes), metabolic process (406 genes), cellular process (502 genes), anatomical structure formation (78 genes) and cellular component organization (134 genes), while the highly expressed genes of female (F2) were mainly enriched on multicellular organismal process (174 genes), localization (147 genes), establishment of localization (130 genes) and biological regulation (250 genes).
KEGG pathway enrichment analysis was also performed in Table S1B. Two-year-old male tongue sole had greater gene expression in genetic information processing, while female’s highly expressed genes were mainly enriched on organismal systems and human diseases pathway.

2.4. Genes with Sex-Biased Expression Potentially Associated with Growth

As a vital mariculture fish, growth differences between male and female tongue sole are attracting more and more attention of researchers. Two-year-old female tongue sole are bigger in size and have higher growth rates than their male counterparts. Strikingly different gene expression patterns were also observed in the brain of the two sexes (Figure 2).
Gene Ontology is known to describe the properties of genes and their products in any organism. Based on GO terms, we attempted to identify sex-biased genes related with sexually dimorphic phenotypic traits. Here, we mainly focus on detecting the sex-biased genes for growth and reproduction associated with GO terms.
According to GO annotation, 5467 genes were annotated by Gene Ontology. In biological process, 115 genes were assigned to growth (GO: 0040007) (Figure 1 and Figure S1). Thirty-three of the genes were sexually differentially expressed, which are called potential growth-associated DEGs in this study (Figure 1 and Table S3).
In two-year-old tongue sole, 30 differentially expressed genes were assigned to growth (Table 2), which accounted for 91% genes of potential growth-associated DEGs. The expression level of these genes in two development stages of tongue sole are plotted in Figure S2, and the significant analysis are listed in Table S3. Twelve genes were significantly highly expressed in M2 compared with F2 (Table 2 and Table S3A), including h3f3b, hsc70, exosc9, vdac2, ppp2r1a, etc. These genes were only highly expressed in M2, with no significant differences in gene expression among M1, F1 and F2 (Table S3A). Another 18 genes were found to be down regulated in two-year-old males compared to the other three samples (Table 2 and Table S3A), including mdk, gap43, ptn, gja1, ninj2, cxcl12b, fgf12, etc. These genes are mainly related to growth, fin regeneration and tissue regeneration. Significantly expressed DEGs were detected in all samples, but interestingly, most of these genes were only up or down regulated in M2 compared to the other samples.

2.5. Genes with Sex-Biased Expression Potentially Associated with Reproduction

Reproduction of tongue sole is another important focus in this study. According to GO annotation, 73 genes were assigned to reproduction (GO: 0000003) in biological process (Figure 1 and Figure S1). Twenty-five of them were differentially expressed between sexes, which were called potential reproduction-associated DEGs in this study (Table 3). The expression level of these genes in two development stages of tongue sole is plotted in Figure S3, and the statistical test analysis is listed in Table S3. For two-year-old fish, 20 genes were found to be up-regulated genes in males, which accounted for the majority of sex-biased genes influencing in reproduction (Table 3 and Table S3B). Those DEGs were only highly expressed in M2, and with no significant differences among the other three samples (Table S3B), which may be closely related to male development process. Nevertheless, only five genes were identified as female-biased genes (Table 3), and these genes were found to be down regulated in M2 compared to F2 (Figure S3B).

2.6. qPCR Validation of RNA-Seq Results

In order to validate the expression pattern of DEGs identified by RNA-Sequencing, we randomly selected eleven genes from DEGs potentially associated with growth or reproduction (Table 2 and Table 3) for qPCR validation, including vdac2, ccnblip1, ropn1l, exosc9, cct7, mdk, ywhaq, gap43, cxcl12b, foxb1 and spin1. M2 was used as the reference in 2−∆∆Ct method to calculate the relative gene expression level. Comparing the relative expression level of eleven selected genes, most results of qPCR were consistent with the results of RNA-Seq (Figure 3). The Pearson’s correlation of log10(fold-change) between qPCR and RNA-Seq was 0.80, indicating the accuracy and reliability of RNA-Seq based transcriptome analysis.

3. Discussion

This study conducted a transcriptome analysis of male and female brain tissues thereby providing the first assessment of the molecular mechanism underlying the reproductive development of tongue sole. We identified sex-biased genes and investigated up/down regulated DEGs with respect to two phenotypic traits (growth and reproduction). Finally, we validated eleven randomly selected potential growth or reproduction associated DEGs by qPCR, thereby supporting the reliability and accuracy of our transcriptome analysis.

3.1. The Brain Transcriptome of Tongue Sole

Brain transcriptome analysis was conducted on male and female tongue sole. The reason we chose one- and two-year-old individuals was two-fold. Firstly, males and females have significant differences in the age of maturity and body weight. Early sexual maturation of males was found at the age of nine months, while 21-month-old females were still immature [25]. The body weight of females was 3.28 times higher than males at 21 months [25]. We tried to capture the transcriptome differences between males and females at a stage with phenotypic divergence. Secondly, we looked for genes that contribute to sexual dimorphism, including upstream genes activity that might contribute to sex determination as well. We chose one- and two-year-old tongue sole for this research. Two-year-old males were mature, while females were still immature. The stage selection was appropriate and effective to characterize the sex-related transcriptomic expression profiles of tongue sole, since we identified a batch of differentially expressed genes in our assembled sequences.

3.2. Male and Female Brain Expression Pattern

Brain plays a central role in coordinating sexual function in mammals, and has remarkable plasticity in teleost [28]. Our results identified more genes to be over expressed in males than females, which is consistent with the findings on Manousaki’s work [20] on sharpsnout seabream, in which the authors identified a seven-fold increase of DEGs in two-year-old fish compared to one-year-old fish. More DEGs were detected in male tongue sole than female during two stages in our study.

3.3. Sexually Differentially Expressed Genes Potentially Associated with Growth

A number of genes involved in growth were found to be down regulated in two-year-old male brains. These included mdk, which encodes a member of a small family of secreted growth factors that binds heparin and promotes cell growth [29,30,31], migration, and angiogenesis [32,33]. It was also reported as a retinoic acid-responsive gene concerned with prenatal development and neuritis growth [34]. In mouse and human embryogenesis, mdk is widely expressed during mid gestation, but expression becomes undetectable in later stages [35,36]. In the brain of adult mice and rats, expression of mdk is reduced to undetectable levels [37]. ptn, often studied together with mdk, is a secreted growth factor that induces neurite outgrowth [38,39,40]. It is also claimed to play essential roles in female reproduction in mice [41]. One of these two genes deficient in mice will result in abnormal reproduction and development [41]. In rodents, these growth factors are highly expressed in early life and decrease to undetectable levels by adulthood in multiple organs [42]. gap43, a growth associated gene [43], plays a key role in nervous system [44] and it is related to nerve development, regeneration, and outgrowth. ninj2 is known to play a role in nerve regeneration and in the formation and function of other tissues [45,46,47]. In this study, the decreased expression level of these genes in M2 might be related to the maturation of male in this stage with growth retardation. However, these genes keep a certain level during two stages in female while it was still growing, which indicated that they might be essential for female growth.
The up regulated genes potentially associated with growth in M2 included h3f3b, which is a replacement histone subtype in active genes [48], involved in DNA damage and cell cycle regulation [49,50]. exosc9 encodes a component of the exosome complex with RNA degradation function [51], which prevents premature differentiation and promotes self-renewal by maintaining a low mRNA level of transcription factor grhl3 [52]. These two genes both have a positive regulatory effect of cell growth. The expression level of h3f3b in two-year-old male (M2) was two times higher than other samples, and exosc9 also shows strikingly higher expression level in two-year-old males (M2), which further suggests an enhanced rate for cell cycle regulation and enhanced self-renewal in the brain of male tongue sole.

3.4. Sexually Differentially Expressed Genes Potentially Associated with Reproduction

A few genes potentially associated with reproduction are discussed as follows. ccnb1ip1, also called as hei10, is a member of the E3 ubiquitin ligase family that acts as a limiting factor for crossing-over during meiosis [53] and involved in male gametogenesis and plays a role in cell cycle regulation [53,54]. When cells enter into mitosis, ccnb1ip1 may interact with Cyclin B1 directly to control protein degradation in human [55], indicating its function of inhibiting cell growth.
ropn1l encodes a sperm protein of the ropporin family, which interacts with Protein Kinase A to regulate glycogen, sugar, and lipid metabolism [56]. Ropporin had been identified as a spermatogenic cell protein, but also has weak expression in normal non-testis tissues in latest studies, including brain tissue [57]. In Tonevitsky’s study [58], ropn1l had an increasing expression level in blood during the relaxation time after 30 min exercise of tested athletes. ropn1l was highly expressed in two-year-old male (M2), indicating its potential function in maturity male tongue sole. cct4, cct5 and cct7 are members of the chaperonin containing TCP1 complex (CCT), act as a molecular chaperone to assist the folding of proteins upon ATP hydrolysis [40,59,60], which are involved in binding of spermatozoa to the zona pellucida. The CCT complex is a large multi-subunit complex that assists the folding of a variety of proteins, including several tubulin and actin proteins [61,62] and it plays a key role in cell cycle regulation, cytoskeleton assembly and cell division [61]. Grantham [61] indicated that full CCT activity is essential for normal cell growth and division in mammalian cell. The functional significance of the increased cct4, cct5 and cct7 expression was unclear, but may have a possible relationship with the quicker maturation in males compared to females.
Other genes found in our research that might have significance but still unknown roles in the brain associated with growth and reproduction. Despite the relatively low number, these genes may become biological significant when interacting with the other genes in functional gene networks controlling and maintaining reproduction and development of tongue sole.

4. Materials and Methods

4.1. Experimental Fish and Sample Collection

All procedures involving sample handling and treatment used in this study were approved by the Animal Care and Use Committee at Heilongjiang River Fisheries Research Institute (ACUC-HRFRI). All experiments were conducted in accordance with the relevant guidelines and regulations. The fish used in this study included 20 male tongue sole (10 one-year-old juveniles and 10 two-year-old adults) and 20 female tongue sole (10 one-year-old juveniles and 10 two-year-old adults). The average body weights of one-year-old males and females were 54 and 99 g, and the average body lengths were 21 and 25 cm, respectively. For two-year-old fish, the average weights were 337 and 1380 g, and average lengths were 27 and 42 cm, for males and females, respectively. Females were much bigger than males. Two-year-old females were still immature, while males were mature. All forty fish came from one full-sib family in National Fish Original Species Farm at Tianjin, China. One-year-old fish were collected on 10 August 2013 and two-year-old on 12 August 2014. The sex of the fish was confirmed anatomically. Tissues of the whole brain from one- and two-year-old fish were collected and stored in 1.5 mL RNAlater tube (QIAGEN, Hilden, Germany), and then stored at −80 °C freezer for subsequent RNA extraction.

4.2. RNA Isolation

Brain samples were taken out of the −80 °C freezer and homogenerated with a TissueRuptor (QIAGEN) to a fine solution and then treated with RNase free DNase I (QIAGEN) for high RNA quality and yield. Total RNA were extracted from each sample using the RNeasy Lipid Tissue Kit (QIAGEN) according to the manufacturer’s instructions. The quality and quantity of RNA were examined by Thermo ScientificTM NanoDropTM 8000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only RNA with OD260/280 ≥ 1.8 and RNA integrity number (RIN) ≥ 7 was stored for the following library construction and sequencing. Equal quantities of high quality RNA from each brain sample were pooled together for subsequent cDNA synthesis and sequencing.

4.3. cDNA Library Construction and Illumina Sequencing

Four RNA-Seq libraries were prepared and sequenced by BerryGenomics sequencing company (Beijing, China). cDNA library was constructed following the protocol of Illumina TrueSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA). The library was sequenced using Illumina HiSeq 2000 platform. In total, 100 bp paired-end reads with an insert size of 330 were obtained. The clean reads from the four transcriptomes were obtained by filtering out low quality reads (reads with quality value of ≤20) and adaptor-only reads from the raw data.

4.4. Sequence Assembly and Analysis

The quality control of RNA-Seq data was conducted by NGS QC Toolkit [63] with default parameters. Clean paired-end reads were aligned to the genome of tongue sole [2] using TopHat [64]. Then the mapping files were analyzed using Cufflinks [65] to assemble the reads into transcripts for each dataset. Complete transcripts were obtained by merging the assemblies of all datasets using Cuffmerge. Duplicate sequences were removed by cd-hit [66].

4.5. Functional Annotation and Ontology

All expressed genes were aligned to NCBI non-redundant (NR) protein database by using BLASTP with E-value < 1 × 10−5. With the result of BLASTP, Gene ontology (GO) [67] annotation was performed using Blast2GO [68]. The top two GO terms associated with transcripts were retrieved, and then the annotation genes were grouped into three categories, biological process, molecular function, and cellular components. Genes annotated on growth (GO: 0040007: BP) and reproduction (GO: 0000003: BP) were selected as potential growth- and reproduction-associated genes, respectively. Pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG) [69] were assigned with differential expressed genes by using the online KEGG Automatic Annotation Server (KAAS) [70].

4.6. Identification of Differentially Expressed Genes (DEGs) and Pathway Enrichment Analysis

FPKM (Fragments Per Kilobase of exon per Million fragments mapped reads) is widely used in RNA-Sequencing to represent the gene expression level [71]. The FPKM is calculated by Cuffdiff [72] method. Then, DEGseq [73] was used to identify sexually differentially expressed genes (DEGs) based on the MARS model (MA-plot-based method with Random Sampling model). Genes with “p-value < 0.001” were selected as DEGs. False discovery rate (FDR) was widely set at 0.05 to determine the threshold for the p-value. Sex-biased genes were strictly identified from DEGs with at least 2-fold difference between genders (|log2(FRKM_ZZ/FPKM_ZW)| ≥ 1). Genes meeting “FPKM_ZZ < 2 and FPKM_ZW < 2” statistical criteria were grouped as non-expressed genes (NEGs). All remaining genes were classified as co-expressed genes (CEGs). This way, all transcribed genes were classified into three groups: DEGs, CEGs and NEGs.
Fisher’s exact test was conducted to further detect functional DEGs for RNA-Seq data [74,75]. The right-tailed Fisher’s exact test was used to identify the enriched GO terms or KEGG categories. All transcribed genes were divided into two sets, differentially expressed genes (DEGs) and non-DEGs. Fisher’s exact test was calculated based on a 2 × 2 contingency table (Table 4). The two columns of the table were separated by DEGs and non-DEGs. The two rows were separated by the given category A and other categories.
The Fisher’s exact test p-value for each category was calculated in the following formula.
p = C a + b a C c + d c C a + b + c + d a + c = ( a + b ) ! ( c + d ) ! ( a + c ) ! ( b + d ) ! ( a + b + c + d ) ! a ! b ! c ! d !
p-Value ≤ 0.05 were selected as a threshold to determine the significant enriched GO terms or KEGG pathways.

4.7. Experimental Validation by qPCR

Sexually differentially expressed genes were validated using quantitative real-time PCR (qPCR). Brain samples were dissected from three male and three female tongue sole at both one and two years of age. Total RNA was extracted from each brain sample and reverse-transcribed using MMLV reverse transcriptase (Invitrogen, Carlsbad, CA, USA). Primers used in qPCR were designed with the help of Primer Premier 6.0. Quantitative real-time PCR was performed using the SYBR Green I Master Mix (TaKaRa, Dalian, China) in Applied Biosystems Prism 7500-fast real-time PCR system. The PCR reaction conditions were as follows: 95 °C for 5 min, followed by 40 amplification cycles of 95 °C for 15 s and 60 °C for 30 s. β-Actin had been validated as a useful internal control for gene expression studies using quantitative real-time PCR in tongue sole [24,76], and was used as the reference gene to normalize the mRNA expression levels of all the validated genes in this study. The sequence of beta-actin gene and primers used in this study were referred from Chen’s study [2]. The fold change of each gene was calculated by the relative quantification (2−∆∆Ct) method [77]. qPCR analysis was repeated three times to confirm gene expression patterns. Expression differences between samples were tested by two-tailed Student’s t-test, and 95% confidence level (p < 0.05) was selected to determine the differentially expressed genes.

5. Conclusions

In this study, we performed a comprehensive analysis of the transcriptome of one- and two-year-old male and female brains of Cynoglossus semilaevis by high-throughput Illumina sequencing. A total of 21,475 unigenes were obtained by transcriptome analysis. Differential expression analysis combined with gene enrichment analysis revealed a number of potential growth and reproduction associated genes. Most of these genes are only up or down regulated in two-year-old males, and had no significant differences in expression among one-year-old males and females and two-year-old females. Potential growth and reproduction associated genes in our work might be the candidate factors affecting sexual dimorphism in tongue sole, laying the foundation to understand the complex process of development and sex determination of this economic valuable species.

Supplementary Materials

Supplementary materials can be found at www.mdpi.com/1422-0067/17/9/1402/s1.

Acknowledgments

We thank Xiaofeng Zhang, Shuqun Xue, and Qiushi Wang for their good suggestion on this manuscript. This research was supported by grants from the Taishan Scholar Leading Talent Team Support Program for Blue Industry of Shandong Province, National Science Foundation of China (31302177), Heilongjiang River Fisheries Research Institute Foundation (HSY201505), and National Science Foundation of Heilongjiang Province (LC201419).

Author Contributions

Xiaowen Sun and Jianguo Lu designed the experiments. Pingping Wang analyzed the data. Min Zheng performed the experiments. Jian Liu and Yongzhuang Liu contributed reagents and materials. Pingping Wang and Jianguo Lu wrote and edited the paper.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

GOGene Ontology
DEGDifferentially Expressed Gene
NEGNon-Expressed Gene
CEGCo-Expressed Gene
FPKMFragments Per Kilobase of exon per Million fragments mapped reads
qPCRquantitative real-time PCR

References

  1. Ma, A.J.; Liu, X.Z.; Xu, Y.J.; Liang, Y.; Zhuang, Z.M. Feeding rhythm and growth of the tongue sole, Cynoglossus semilaevis Gunther, during its early life stages. Aquac. Res. 2006, 37, 586–593. [Google Scholar]
  2. Chen, S.; Zhang, G.; Shao, C.; Huang, Q.; Liu, G.; Zhang, P.; Song, W.; An, N.; Chalopin, D.; Volff, J.N.; et al. Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. Nat. Genet. 2014, 46, 253–260. [Google Scholar] [CrossRef] [PubMed]
  3. Sun, Y.Y.; Zhang, Q.Q.; Qi, J.; Chen, Y.J.; Zhong, Q.W.; Li, C.M.; Yu, Y.; Li, S.; Wang, Z.G. Identification of differential genes in the ovary relative to the testis and their expression patterns in half-smooth tongue sole (Cynoglossus semilaievis). J. Genet. Genom. 2010, 37, 137–145. [Google Scholar] [CrossRef]
  4. Baroiller, J.F.; D’Cotta, H.; Saillant, E. Environmental effects on fish sex determination and differentiation. Sex. Dev. 2009, 3, 118–135. [Google Scholar] [CrossRef] [PubMed]
  5. Ospina-Alvarez, N.; Piferrer, F. Temperature-dependent sex determination in fish revisited: Prevalence, a single sex ratio response pattern, and possible effects of climate change. PLoS ONE 2008, 3, e2837. [Google Scholar] [CrossRef] [PubMed]
  6. Shang, E.H.H.; Yu, R.M.K.; Wu, R.S.S. Hypoxia affects sex differentiation and development, leading to a male-dominated population in zebrafish (Danio rerio). Environ. Sci. Technol. 2006, 40, 3118–3122. [Google Scholar] [PubMed]
  7. Francis, R.C. Sexual lability in teleosts: Developmental factors. Q. Rev. Biol. 1992, 67, 18. [Google Scholar] [CrossRef]
  8. Weltzien, F.A.; Andersson, E.; Andersen, O.; Shalchian-Tabrizi, K.; Norberg, B. The brain-pituitary-gonad axis in male teleosts, with special emphasis on flatfish (Pleuronectiformes). Comp. Biochem. Physiol. Part A Mol. Integr. Physiol. 2004, 137, 447–477. [Google Scholar] [CrossRef] [PubMed]
  9. Gahr, M. Male Japanese quails with female brains do not show male sexual behaviors. Proc. Natl. Acad. Sci. USA 2003, 100, 7959–7964. [Google Scholar] [CrossRef] [PubMed]
  10. Dennis, C. Brain development: The most important sexual organ. Nature 2004, 427, 390–392. [Google Scholar] [CrossRef]
  11. Davies, W.; Wilkinson, L.S. It is not all hormones: Alternative explanations for sexual differentiation of the brain. Brain Res. 2006, 1126, 36–45. [Google Scholar] [CrossRef] [PubMed]
  12. Dewing, P.; Shi, T.; Horvath, S.; Vilain, E. Sexually dimorphic gene expression in mouse brain precedes gonadal differentiation. Brain Res. Mol. Brain Res. 2003, 118, 82–90. [Google Scholar] [CrossRef]
  13. Yang, X.; Schadt, E.E.; Wang, S.; Wang, H.; Arnold, A.P.; Ingram-Drake, L.; Drake, T.A.; Lusis, A.J. Tissue-specific expression and regulation of sexually dimorphic genes in mice. Genome Res. 2006, 16, 995–1004. [Google Scholar] [CrossRef] [PubMed]
  14. Parsch, J.; Ellegren, H. The evolutionary causes and consequences of sex-biased gene expression. Nat. Rev. Gene. 2013, 14, 83–87. [Google Scholar] [CrossRef] [PubMed]
  15. Khatri, P.; Draghici, S.; Ostermeier, G.C.; Krawetz, S.A. Profiling gene expression using onto-express. Genomics 2002, 79, 266–270. [Google Scholar] [CrossRef] [PubMed]
  16. Parisi, M.; Nuttall, R.; Edwards, P.; Minor, J.; Naiman, D.; Lu, J.; Doctolero, M.; Vainer, M.; Chan, C.; Malley, J.; et al. A survey of ovary-, testis-, and soma-biased gene expression in Drosophila melanogaster adults. Genome Biol. 2004, 5, R40. [Google Scholar] [CrossRef] [PubMed][Green Version]
  17. Small, C.L.; Shima, J.E.; Uzumcu, M.; Skinner, M.K.; Griswold, M.D. Profiling gene expression during the differentiation and development of the murine embryonic gonad. Biol. Reprod. 2005, 72, 492–501. [Google Scholar] [CrossRef] [PubMed]
  18. Santos, E.M.; Workman, V.L.; Paull, G.C.; Filby, A.L.; van Look, K.J.; Kille, P.; Tyler, C.R. Molecular basis of sex and reproductive status in breeding zebrafish. Physiol. Genom. 2007, 30, 111–122. [Google Scholar] [CrossRef] [PubMed]
  19. Tao, W.; Yuan, J.; Zhou, L.; Sun, L.; Sun, Y.; Yang, S.; Li, M.; Zeng, S.; Huang, B.; Wang, D. Characterization of gonadal transcriptomes from Nile tilapia (Oreochromis niloticus) reveals differentially expressed genes. PLoS ONE 2013, 8, e63604. [Google Scholar] [CrossRef] [PubMed]
  20. Manousaki, T.; Tsakogiannis, A.; Lagnel, J.; Sarropoulou, E.; Xiang, J.Z.; Papandroulakis, N.; Mylonas, C.C.; Tsigenopoulos, C.S. The sex-specific transcriptome of the hermaphrodite sparid sharpsnout seabream (Diplodus puntazzo). BMC Genom. 2014, 15, 1–16. [Google Scholar] [CrossRef] [PubMed]
  21. Sharma, E.; Kunstner, A.; Fraser, B.A.; Zipprich, G.; Kottler, V.A.; Henz, S.R.; Weigel, D.; Dreyer, C. Transcriptome assemblies for studying sex-biased gene expression in the guppy, Poecilia reticulata. BMC Genom. 2014, 15, 400. [Google Scholar] [CrossRef] [PubMed]
  22. Lu, J.; Luan, P.; Zhang, X.; Xue, S.; Peng, L.; Mahbooband, S.; Sun, X. Gonadal transcriptomic analysis of yellow catfish (Pelteobagrus fulvidraco): Identification of sex-related genes and genetic markers. Physiol. Genom. 2014, 46, 798–807. [Google Scholar] [CrossRef] [PubMed]
  23. Ji, X.S.; Chen, S.L.; Jiang, Y.L.; Xu, T.J.; Yang, J.F.; Tian, Y.S. Growth differences and differential expression analysis of pituitary adenylate cyclase activating polypeptide (PACAP) and growth hormone-releasing hormone (GHRH) between the sexes in half-smooth tongue sole Cynoglossus semilaevis. Gen. Comp. Endocrinol. 2011, 170, 99–109. [Google Scholar] [CrossRef] [PubMed]
  24. Ma, Q.; Liu, S.F.; Zhuang, Z.M.; Lin, L.; Sun, Z.Z.; Liu, C.L.; Su, Y.Q.; Tang, Q.S. Genomic structure, polymorphism and expression analysis of growth hormone-releasing hormone and pituitary adenylate cyclase activating polypeptide genes in the half-smooth tongue sole (Cynoglossus semilaevis). Genet. Mol. Res. 2011, 10, 3828–3846. [Google Scholar] [CrossRef] [PubMed]
  25. Ji, X.S.; Liu, H.W.; Chen, S.L.; Jiang, Y.L.; Tian, Y.S. Growth differences and dimorphic expression of growth hormone (GH) in female and male Cynoglossus semilaevis after male sexual maturation. Mar. Genom. 2011, 4, 9–16. [Google Scholar] [CrossRef] [PubMed]
  26. Ma, Q.; Liu, S.F.; Zhuang, Z.M.; Lin, L.; Sun, Z.Z.; Liu, C.L.; Ma, H.; Su, Y.Q.; Tang, Q.S. Genomic structure, polymorphism and expression analysis of the growth hormone (GH) gene in female and male Half-smooth tongue sole (Cynoglossus semilaevis). Gene 2012, 493, 92–104. [Google Scholar] [CrossRef] [PubMed]
  27. Wang, Z.; Gerstein, M.; Snyder, M. RNA-Seq: A revolutionary tool for transcriptomics. Nat. Rev. Genet. 2009, 10, 57–63. [Google Scholar] [CrossRef] [PubMed]
  28. Trabzuni, D.; Ramasamy, A.; Imran, S.; Walker, R.; Smith, C.; Weale, M.E.; Hardy, J.; Ryten, M.; North American Brain Expression Consortium. Widespread sex differences in gene expression and splicing in the adult human brain. Nat. Commun. 2013, 4, 2771. [Google Scholar] [CrossRef] [PubMed]
  29. Sakakima, H.; Yoshida, Y.; Yamazaki, Y.; Matsuda, F.; Ikutomo, M.; Ijiri, K.; Muramatsu, H.; Muramatsu, T.; Kadomatsu, K. Disruption of the Midkine gene (MDK) delays degeneration and regeneration in injured peripheral nerve. J. Neurosci. Res. 2009, 87, 2908–2915. [Google Scholar] [CrossRef] [PubMed]
  30. Singec, I.; Crain, A.M.; Tobe, B.T.D.; Hou, J.; Talantova, M.; Doctor, K.S.; Choi, J.; Huang, X.; Gutierrez, G.J.; Wolf, D.A.; et al. Midkine (MDK): A newly-recognized endogenously-produced cytokine that, through an autocrine mechanism, is pivotal for neural induction. Neuroreport 2014, 25, 150. [Google Scholar]
  31. Taylor, T.D.; Noguchi, H.; Totoki, Y.; Toyoda, A.; Kuroki, Y.; Dewar, K.; Lloyd, C.; Itoh, T.; Takeda, T.; Kim, D.W.; et al. Human chromosome 11 DNA sequence and analysis including novel gene identification. Nature 2006, 440, 497–500. [Google Scholar] [CrossRef] [PubMed]
  32. Muramatsu, T. Midkine: A Promising Molecule for Drug Development to Treat Diseases of the Central Nervous System. Curr. Pharm. Des. 2011, 17, 410–423. [Google Scholar] [CrossRef] [PubMed]
  33. Weckbach, L.T.; Muramatsu, T.; Walzog, B. Midkine in Inflammation. Sci. World J. 2011, 11, 2491–2505. [Google Scholar] [CrossRef] [PubMed]
  34. Kaname, T.; Kuwano, A.; Murano, I.; Uehara, K.; Muramatsu, T.; Kajii, T. Midkine gene (MDK), a gene for prenatal differentiation and neuroregulation, maps to band 11p11.2 by fluorescence in situ hybridization. Genomics 1993, 17, 514–515. [Google Scholar] [CrossRef] [PubMed]
  35. Kadomatsu, K.; Huang, R.P.; Suganuma, T.; Murata, F.; Muramatsu, T. A retinoic acid responsive gene MK found in the teratocarcinoma system is expressed in spatially and temporally controlled manner during mouse embryogenesis. J. Cell Biol. 1990, 110, 607–616. [Google Scholar] [CrossRef] [PubMed]
  36. Fan, Q.W.; Muramatsu, T.; Kadomatsu, K. Distinct expression of midkine and pleiotrophin in the spinal cord and placental tissues during early mouse development. Dev. Growth Differ. 2000, 42, 113–119. [Google Scholar] [CrossRef] [PubMed]
  37. Kretschmer, P.J.; Fairhurst, J.L.; Decker, M.M.; Chan, C.P.; Gluzman, Y.; Bohlen, P.; Kovesdi, I. Cloning, characterization and developmental regulation of two members of a novel human gene family of neurite outgrowth-promoting proteins. Growth Factors 1991, 5, 99–114. [Google Scholar] [CrossRef] [PubMed]
  38. Sortwell, C.E.; Collier, T.J.; Terpstra, B.T.; Thompson, V.B.; O’Malley, J.; Steece-Collier, K.; Wohlgenant, S.M.; Daley, B.F.; Mandel, R.J.; Lipton, J.W. The potential for the trophic factor pleiotrophin (PTN) to protect the Parkinsonian brain. Cell Transplant. 2008, 17, 481. [Google Scholar]
  39. Xu, C.Y.; Zhu, S.Y.; Wu, M.Y.; Han, W.; Yu, Y. Functional receptors and intracellular signal pathways of midkine (MK) and pleiotrophin (PTN). Biol. Pharm. Bull. 2014, 37, 511–520. [Google Scholar] [CrossRef] [PubMed]
  40. Human Genome Sequencing, C. Finishing the euchromatic sequence of the human genome. Nature 2004, 431, 931–945. [Google Scholar]
  41. Muramatsu, H.; Zou, P.; Kurosawa, N.; Ichihara-Tanaka, K.; Maruyama, K.; Inoh, K.; Sakai, T.; Chen, L.; Sato, M.; Muramatsu, T. Female infertility in mice deficient in midkine and pleiotrophin, which form a distinct family of growth factors. Genes Cells 2006, 11, 1405–1417. [Google Scholar] [CrossRef] [PubMed]
  42. Jee, Y.H.; Lebenthal, Y.; Chaemsaithong, P.; Yan, G.; Peran, I.; Wellstein, A.; Romero, R.; Baron, J. Midkine and pleiotrophin concentrations in amniotic fluid in healthy and complicated pregnancies. PLoS ONE 2016, 11, e0153325. [Google Scholar] [CrossRef] [PubMed]
  43. De Groen, P.C.; Eggen, B.J.; Gispen, W.H.; Schotman, P.; Schrama, L.H. Cloning and promoter analysis of the human B-50/GAP-43 gene. J. Mol. Neurosci. 1995, 6, 109–119. [Google Scholar] [CrossRef] [PubMed]
  44. Benowitz, L.I.; Routtenberg, A. GAP-43: An intrinsic determinant of neuronal development and plasticity. Trends Neurosci. 1997, 20, 84–91. [Google Scholar] [CrossRef]
  45. Scherer, S.E.; Muzny, D.M.; Buhay, C.J.; Chen, R.; Cree, A.; Ding, Y.; Dugan-Rocha, S.; Gill, R.; Gunaratne, P.; Harris, R.A.; et al. The finished DNA sequence of human chromosome 12. Nature 2006, 440, 346–351. [Google Scholar] [CrossRef] [PubMed]
  46. Araki, T.; Milbrandt, J. Ninjurin2, a novel homophilic adhesion molecule, is expressed in mature sensory and enteric neurons and promotes neurite outgrowth. J. Neurosci. 2000, 20, 187–195. [Google Scholar] [PubMed]
  47. Bis, J.C.; DeStefano, A.; Liu, X.M.; Brody, J.A.; Choi, S.H.; Verhaaren, B.F.J.; Debette, S.; Ikram, M.A.; Shahar, E.; Butler, K.R.; et al. Associations of NINJ2 sequence variants with incident ischemic stroke in the cohorts for heart and aging in genomic epidemiology (charge) consortium. PLoS ONE 2014, 9, e99798. [Google Scholar] [CrossRef] [PubMed][Green Version]
  48. Albig, W.; Bramlage, B.; Gruber, K.; Klobeck, H.G.; Kunz, J.; Doenecke, D. The human replacement histone H3.3B gene (H3F3B). Genomics 1995, 30, 264–272. [Google Scholar] [CrossRef] [PubMed]
  49. De Luca, P.; Vazquez, E.S.; Moiola, C.P.; Zalazar, F.; Cotignola, J.; Gueron, G.; Gardner, K.; de Siervi, A. BRCA1 loss induces GADD153-mediated doxorubicin resistance in prostate cancer. Mol. Cancer Res. 2011, 9, 1078–1090. [Google Scholar] [CrossRef] [PubMed]
  50. De Luca, P.; Moiola, C.P.; Zalazar, F.; Gardner, K.; Vazquez, E.S.; de Siervi, A. BRCA1 and p53 regulate critical prostate cancer pathways. Prostate Cancer Prostatic Dis. 2013, 16, 233–238. [Google Scholar] [CrossRef] [PubMed]
  51. Van Dijk, E.L.; Schilders, G.; Pruijn, G.J. Human cell growth requires a functional cytoplasmic exosome, which is involved in various mRNA decay pathways. RNA 2007, 13, 1027–1035. [Google Scholar] [CrossRef] [PubMed]
  52. Mistry, D.S.; Chen, Y.F.; Sen, G.L. Progenitor function in self-renewing human epidermis is maintained by the exosome. Cell Stem Cell 2012, 11, 127–135. [Google Scholar] [CrossRef] [PubMed]
  53. Strong, E.R.; Schimenti, J.C. Evidence implicating CCNB1IP1, a RING domain-containing protein required for meiotic crossing over in mice, as an E3 SUMO ligase. Genes 2010, 1, 440–451. [Google Scholar] [CrossRef] [PubMed]
  54. Gronholm, M.; Muranen, T.; Toby, G.G.; Utermark, T.; Hanemann, C.O.; Golemis, E.A.; Carpen, O. A functional association between merlin and HEI10, a cell cycle regulator. Oncogene 2006, 25, 4389–4398. [Google Scholar] [CrossRef] [PubMed]
  55. Ma, Z.; Lin, M.; Li, K.; Fu, Y.; Liu, X.; Yang, D.; Zhao, Y.; Zheng, J.; Sun, B. Knocking down SMC1A inhibits growth and leads to G2/M arrest in human glioma cells. Int. J. Clin. Exp. Pathol. 2013, 6, 862–869. [Google Scholar] [PubMed]
  56. Chen, L.; Kass, R.S. A-kinase anchoring proteins: Different partners, different dance. Nat. Cell Biol. 2005, 7, 1050–1051. [Google Scholar] [CrossRef] [PubMed]
  57. Chen, J.; Cai, Z.M.; Gui, Y.T. Advances in the researches of spermatogenic protein, Ropporin. Natl. J. Androl. 2009, 15, 3. [Google Scholar]
  58. Tonevitsky, A.G.; Maltseva, D.V.; Abbasi, A.; Samatov, T.R.; Sakharov, D.A.; Shkurnikov, M.U.; Lebedev, A.E.; Galatenko, V.V.; Grigoriev, A.I.; Northoff, H. Dynamically regulated miRNA-mRNA networks revealed by exercise. BMC Physiol. 2013, 13, 11. [Google Scholar] [CrossRef] [PubMed]
  59. Sergeeva, O.A.; Chen, B.; Haase-Pettingell, C.; Ludtke, S.J.; Chiu, W.; King, J.A. Human CCT4 and CCT5 chaperonin subunits expressed in Escherichia coli form biologically active homo-oligomers. J. Biol. Chem. 2013, 288, 17734–17744. [Google Scholar] [CrossRef] [PubMed]
  60. Damours, O.; Sullivan, R. Evaluating Fertility of Mammalian Spermatozoa e.g., Bovine Spermatozoa Involves Assessing in Spermatozoa Sample such as Bovine Sperm, Amount or Activity of Specific Sperm Fertility Protein, Where Amount/Activity Is Predictive of Fertility. WO2010025548-A1, 11 March 2010. [Google Scholar]
  61. Grantham, J.; Brackley, K.I.; Willison, K.R. Substantial CCT activity is required for cell cycle progression and cytoskeletal organization in mammalian cells. Exp. Cell Res. 2006, 312, 2309–2324. [Google Scholar] [CrossRef] [PubMed]
  62. Feldman, D.E.; Thulasiraman, V.; Ferreyra, R.G.; Frydman, J. Formation of the VHL-elongin BC tumor suppressor complex is mediated by the chaperonin TRiC. Mol. Cell 1999, 4, 1051–1061. [Google Scholar] [CrossRef]
  63. Patel, R.K.; Jain, M. NGS QC Toolkit: A toolkit for quality control of next generation sequencing data. PLoS ONE 2012, 7, e30619. [Google Scholar] [CrossRef] [PubMed]
  64. Kim, D.; Salzberg, S.L. TopHat-Fusion: An algorithm for discovery of novel fusion transcripts. Genome Biol. 2011, 12. [Google Scholar] [CrossRef] [PubMed]
  65. Roberts, A.; Pimentel, H.; Trapnell, C.; Pachter, L. Identification of novel transcripts in annotated genomes using RNA-Seq. Bioinformatics 2011, 27, 2325–2329. [Google Scholar] [CrossRef] [PubMed]
  66. Fu, L.; Niu, B.; Zhu, Z.; Wu, S.; Li, W. CD-HIT: Accelerated for clustering the next-generation sequencing data. Bioinformatics 2012, 28, 3150–3152. [Google Scholar] [CrossRef] [PubMed]
  67. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene Ontology: Tool for the unification of biology. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed]
  68. Conesa, A.; Gotz, S.; Garcia-Gomez, J.M.; Terol, J.; Talon, M.; Robles, M. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21, 3674–3676. [Google Scholar] [CrossRef] [PubMed]
  69. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef] [PubMed]
  70. Moriya, Y.; Itoh, M.; Okuda, S.; Yoshizawa, A.C.; Kanehisa, M. KAAS: An automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007, 35, 182–185. [Google Scholar] [CrossRef] [PubMed]
  71. Mortazavi, A.; Williams, B.A.; Mccue, K.; Schaeffer, L.; Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 2008, 5, 621–628. [Google Scholar] [CrossRef] [PubMed]
  72. Trapnell, C.; Hendrickson, D.G.; Sauvageau, M.; Goff, L.; Rinn, J.L.; Pachter, L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat. Biotechnol. 2013, 31, 46–53. [Google Scholar] [CrossRef] [PubMed]
  73. Wang, L.; Feng, Z.; Wang, X.; Wang, X.; Zhang, X. DEGseq: An R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 2010, 26, 136–138. [Google Scholar] [CrossRef] [PubMed]
  74. Robinson, M.D.; Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11, R25. [Google Scholar] [CrossRef] [PubMed]
  75. Bloom, J.S.; Khan, Z.; Kruglyak, L.; Singh, M.; Caudy, A.A. Measuring differential gene expression by short read sequencing: Quantitative comparison to 2-channel gene expression microarrays. BMC Genom. 2009, 10, 221. [Google Scholar] [CrossRef] [PubMed]
  76. Li, Z.; Yang, L.; Wang, J.; Shi, W.; Pawar, R.A.; Liu, Y.; Xu, C.; Cong, W.; Hu, Q.; Lu, T.; et al. β-Actin is a useful internal control for tissue-specific gene expression studies using quantitative real-time PCR in the half-smooth tongue sole Cynoglossus semilaevis challenged with LPS or Vibrio anguillarum. Fish Shellfish Immunol. 2010, 29, 89–93. [Google Scholar] [CrossRef] [PubMed]
  77. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Work-flow of the transcriptome analysis of tongue sole.
Figure 1. Work-flow of the transcriptome analysis of tongue sole.
Ijms 17 01402 g001
Figure 2. Gene expression profiling in brain tissue between male and female tongue sole at two developmental stages. Non-expressed genes (NEGs) marked as black dots. Differentially expressed genes (DEGs) marked as dark blue and light blue dots. Co-expressed genes (CEGs) marked as red dots.
Figure 2. Gene expression profiling in brain tissue between male and female tongue sole at two developmental stages. Non-expressed genes (NEGs) marked as black dots. Differentially expressed genes (DEGs) marked as dark blue and light blue dots. Co-expressed genes (CEGs) marked as red dots.
Ijms 17 01402 g002
Figure 3. Brain transcriptome validation by qPCR using eleven selected genes: (A) Gene relative expression in qPCR where Error bars represent standard deviation of technical replicates and the expression differences between M2 and F2 were tested by two-tailed Student’s t-test, with 95% confidence level (p < 0.05) being labeled with an asterisk; and (B) gene expression in RNA-Seq.
Figure 3. Brain transcriptome validation by qPCR using eleven selected genes: (A) Gene relative expression in qPCR where Error bars represent standard deviation of technical replicates and the expression differences between M2 and F2 were tested by two-tailed Student’s t-test, with 95% confidence level (p < 0.05) being labeled with an asterisk; and (B) gene expression in RNA-Seq.
Ijms 17 01402 g003
Table 1. The Gene ontology enrichment analysis of differential expressed genes. For Fisher’s exact test, p-value ≤ 0.05 is labeled as yellow, and p-value ≤ 0.01 is labeled as red.
Table 1. The Gene ontology enrichment analysis of differential expressed genes. For Fisher’s exact test, p-value ≤ 0.05 is labeled as yellow, and p-value ≤ 0.01 is labeled as red.
GO CategoryGO TermTotal No. of GenesM1 vs. F1M2 vs. F2
M1F1M2F2
Cellular ComponentGO:0005576 extracellular region3101094049
GO:0005623 cell32386938425400
GO:0019012 virion113052
GO:0031974 membrane-enclosed lumen408828331
GO:0031975 envelope105702218
GO:0032991 macromolecular complex10143426175151
GO:0043226 organelle16693528282212
GO:0044421 extracellular region part214343433
GO:0044422 organelle part10082915193127
GO:0044423 virion part113052
GO:0044456 synapse part7531323
GO:0044464 cell part32386938425400
GO:0045202 synapse10261528
Molecular FunctionGO:0003824 catalytic activity21053413298210
GO:0005198 structural molecule activity1129182042
GO:0005215 transporter activity3931953068
GO:0005488 binding27805331380318
GO:0009055 electron carrier activity60011
GO:0015457 auxiliary transport protein activity70112
GO:0016209 antioxidant activity140012
GO:0030234 enzyme regulator activity200042235
GO:0030528 transcription regulator activity188211414
GO:0042056 chemoattractant activity20001
GO:0045182 translation regulator activity261082
GO:0045499 chemorepellent activity10000
GO:0060089 molecular transducer activity430321935
Biological ProcessGO:0000003 reproduction7300205
GO:0001906 cell killing30001
GO:0002376 immune system process189212826
GO:0008152 metabolic process27564637406297
GO:0008371 obsolete biological process584048
GO:0009987 cellular process35537148502403
GO:0010926 anatomical structure formation378887843
GO:0016032 viral reproduction301188
GO:0016043 cellular component organization65215813477
GO:0016265 death169223214
GO:0022414 reproductive process72111910
GO:0022610 biological adhesion113101315
GO:0032501 multicellular organismal process13153417151174
GO:0032502 developmental process11242614147129
GO:0040007 growth115531218
GO:0040011 locomotion152111822
GO:0043473 pigmentation17843326198216
GO:0044085 cellular component biogenesis32711227558
GO:0048511 rhythmic process110001
GO:0050896 response to stimulus6061998183
GO:0051179 localization9553011120147
GO:0051234 establishment of localization7903010100130
GO:0051704 multi-organism process55211111
GO:0065007 biological regulation19273829219250
All genes annotated by GO546710463711619
Table 2. Differentially expressed genes (DEGs) in two-year-old tongue sole potentially associated with growth.
Table 2. Differentially expressed genes (DEGs) in two-year-old tongue sole potentially associated with growth.
GeneFPKMp-ValueAnnotationGO Function
M2F2M2 vs. F2
DEGs Up Regulated in M2
h3f3b670.6193.21.9E-56Histone H3.3Positive regulation of cell growth (GO:0030307)
hsc70247.90.15.5E-34Heat shock cognate 70 kDa proteinFin regeneration (GO:0031101)
exosc9164.613.21.7E-32Exosome complex component RRP45Positive regulation of cell growth (GO:0030307)
vdac2164.451.17.9E-14Voltage-dependent anion-selective channel protein 2Fin regeneration (GO:0031101)
ppp2r1a237.7108.19.4E-11Serine/threonine-protein phosphatase 2A 65 kDa regulatory subunit A β isoformGrowth (GO:0040007)
timm5044.612.12.5E-05Mitochondrial import inner membrane translocase subunit TIM50Growth (GO:0040007)
bmp230.66.15.3E-05Bone morphogenetic protein 2Growth (GO:0040007)
ascl117.01.19.1E-05Achaete-scute homolog 1Sensory epithelium regeneration (GO:0070654)
suv420h131.87.31.1E-04Histone-lysine N-methyltransferase SUV420H1Regulation of multicellular organism growth (GO:0040014)
orc115.71.11.8E-04Origin recognition complex subunit 1Regulation of multicellular organism growth (GO:0040014)
nenf26.85.72.2E-04Neudesin neurotrophic factorGrowth (GO:0040007)
dnajc242.616.09.6E-04dnaJ homolog subfamily C member 2Negative regulation of cell growth (GO:0030308)
DEGs Down Regulated in M2
mdk5.3667.81.7E-157MidkineGrowth (GO:0040007)
gap432.3442.42.3E-100Growth associated protein 43Tissue regeneration (GO:0042246)
ptn30.3287.13.1E-56PleiotrophinGrowth (GO:0040007)
gja11.9222.06.1E-54Gap junction α-1 proteinFin regeneration (GO:0031101)
ninj255.5255.51.1E-34Ninjurin-2Tissue regeneration (GO:0042246)
cxcl12b33.3131.81.1E-16Stromal cell-derived factor 1 precursorFin regeneration (GO:0031101)
fgf122.056.15.0E-15Fibroblast growth factor 12Growth (GO:0040007)
gcfc223.4102.52.3E-14GC-rich sequence DNA-binding factor 2Fin regeneration (GO:0031101)
sox20.738.55.3E-11Transcription factor SOX-2Fin regeneration (GO:0031101)
fyna2.838.21.0E-09Tyrosine-protein kinase fynaFin regeneration (GO:0031101)
epha41.734.51.9E-09Ephrin type-A receptor 4Negative regulation of collateral sprouting of intact axon in response to injury (GO:0048685)
bdnf1.929.36.6E-08Brain-derived neurotrophic factorGrowth (GO:0040007)
cdh20.624.81.5E-07cadherin-2-likeTissue regeneration (GO:0042246)
rerg3.029.92.6E-07Ras-related and estrogen-regulated growth inhibitorNegative regulation of cell growth (GO:0030308)
ctgf2.522.99.4E-06Connective tissue growth factorRegulation of cell growth (GO:0001558)
igfbp10.217.99.8E-06Insulin-like growth factor-binding protein 1Regulation of cell growth (GO:0001558)
klf62.622.51.3E-05Krueppel-like factor 6Organ growth (GO:0035265)
sptbn17.027.51.7E-04Spectrin β chain, non-erythrocytic 1Positive regulation of multicellular organism growth (GO:0040018)
Table 3. Differentially expressed genes (DEGs) in two-year-old tongue sole potentially associated with reproduction.
Table 3. Differentially expressed genes (DEGs) in two-year-old tongue sole potentially associated with reproduction.
GeneFPKMp-ValueAnnotationGO Function
M2F2M2 vs. F2
pleaDEGs Up Regulated in M2
ccnb1ip1661.42.73.6E-138E3 ubiquitin-protein ligase CCNB1IP1Reproductive cellular process (GO:0048610)
rsph1516.75.71.4E-118Radial spoke head 1 homologReproduction (GO:0000003)
ropn1l343.40.81.5E-66Ropporin-1-like proteinSpermatid development (GO:0007286)
alkbh5180.517.94.5E-33RNA demethylase ALKBH5Spermatogenesis (GO:0007283)
henmt1150.30.02.4E-29Small RNA 2’-O-methyltransferaseOocyte development (GO:0048599)
cct5358.3136.83.4E-21Chaperonin containing T-complex polypeptide 5Binding of sperm to zona pellucida (GO:0007339)
ovol1112.10.13.2E-19Putative transcription factor Ovo-like 1Spermatogenesis (GO:0007283)
klhl1054.70.07.0E-14Kelch-like protein 10Spermatid development (GO:0007286)
cct7329.3153.71.4E-13Cct7 proteinBinding of sperm to zona pellucida (GO:0007339)
cct4305.8147.98.8E-12Cct4 proteinBinding of sperm to zona pellucida (GO:0007339)
fem1b44.63.02.5E-10Protein fem-1 homolog B-likeEpithelial cell maturation involved in prostate gland development (GO:0060743)
sf1103.433.79.2E-09Splicing factor 1Leydig cell differentiation (GO:0033327)
styx69.820.23.6E-07Serine/threonine/tyrosine-interacting proteinSpermatogenesis (GO:0007283)
ptbp145.39.41.2E-06Polypyrimidine tract-binding protein 1-likePositive regulation of cortical granule exocytosis by elevation of cytosolic calcium ion concentration (GO:0060472)
tbp41.89.04.4E-06TATA box binding proteinSpermatogenesis (GO:0007283)
pld624.92.61.1E-05Phospholipase D family member 6P granule organization (GO:0030719)
polr2h50.915.83.2E-05DNA-directed RNA polymerases I%2C II%2C and III subunit RPABC3Positive regulation of viral transcription (GO:0050434)
pacrg38.810.57.8E-05Parkin coregulated gene proteinSpermatid development (GO:0007286)
tceb1137.675.31.2E-04Transcription elongation factor B polypeptide 1Positive regulation of viral transcription (GO:0050434)
tdrd913.30.62.9E-04Putative ATP-dependent RNA helicase TDRD9Fertilization (GO:0009566); DNA methylation during gametogenesis (GO:0043046)
DEGs Down Regulated in M2
ywhaq63.5229.41.1E-25YWHAQ proteinGermarium-derived oocyte fate determination (GO:0007294)
cxcl12b33.3131.81.1E-16Stromal cell-derived factor 1-likeGerm cell migration (GO:0008354)
spin112.444.64.0E-06Spindlin-1-likeGamete generation (GO:0007276)
stau241.780.21.1E-04Double-stranded RNA-binding protein Staufen homolog 2Germ cell migration (GO:0008354)
sptbn17.027.51.7E-04Spectrin β chain, non-erythrocytic 1Fertilization (GO:0009566)
Table 4. 2 × 2 contingency table for Fisher’s exact test.
Table 4. 2 × 2 contingency table for Fisher’s exact test.
Gene AnnotationDEGsNon-DEGsAll Transcribed Genes
Category Aaba + b
Other categoriescdc + d
All categoriesa + cb + da + b + c + d
“a” represented the No. of DEGs annotated on category A; “b” represented the No. of non-DEGs annotated on category A; “c” represented the No. of DEGs annotated on the other categories except category A; and “d” represented the No. of non-DEGs annotated on the other categories except category A.
Back to TopTop