Transcriptomic Profiles of Brain Provide Insights into Molecular Mechanism of Feed Conversion Efficiency in Crucian Carp (Carassius auratus)

Feed efficiency is an economically crucial trait for cultured animals, however, progress has been scarcely made in the genetic analyses of feed conversion efficiency (FCE) in fish because of the difficulties in measurement of trait phenotypes. In the present investigation, we present the first application of RNA sequencing (RNA-Seq) combined with differentially expressed genes (DEGs) analysis for identification of functional determinants related to FCE at the gene level in an aquaculture fish, crucian carp (Carassius auratus). Brain tissues of six crucian carp with extreme FCE performances were subjected to transcriptome analysis. A total of 544,612 unigenes with a mean size of 644.38 bp were obtained from Low- and High-FCE groups, and 246 DEGs that may be involved in FCE traits were identified in these two groups. qPCR confirmed that genes previously identified as up- or down-regulated by RNA-Seq were effectively up- or down-regulated under the studied conditions. Thirteen key genes, whose functions are associated with metabolism (Dgkk, Mgst3 and Guk1b), signal transduction (Vdnccsa1b, Tgfα, Nr4a1 and Tacr2) and growth (Endog, Crebrtc2, Myh7, Myh1, Myh14 and Igfbp7) were identified according to GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) annotations. Our novel findings provide useful pathway information and candidate genes for future studies of genetic mechanisms underlying FCE in crucian carp.


Introduction
Aquaculture products are among the most important sources of nutrition for people around the world [1]. Fish is a key source of protein, essential amino-acids and minerals, especially in low-income, food-deficient countries [2][3][4]. In the past decades, yield-related traits in fish have been extensively studied by researchers to meet the ever-increasing global demand and increase financial returns [5,6].
Since feed cost comprises about 65-75% of the total production cost in most aquaculture species [7], the increasing production demand with huge input cost kept up pressures on us to investigate how to improve feed efficiency. In this sense, breeding fish with high feed efficiency would mean big savings and enhance the profitability of producers.
Feed efficiency was usually measured as feed conversion ratio (FCR), which is the ratio of feed intake to body weight gain [8]. To date, some DNA variants playing roles in feed efficiency have been proposed in livestock [8][9][10][11] and poultry industry [12][13][14][15] using quantitative trait loci (QTL) mapping and association studies. Those identified candidate genes pertain to numerous biological processes [10,16], suggesting that biological strategies involving feed efficiency are diverse.
RNA sequencing (RNA-Seq) technology combined with the analysis of differentially expressed genes (DEGs) are reliable and precise ways to provide useful information for molecular mechanism of complex characters [17]. Moreover, RNA-Seq has attracted considerable interests and received great success concerning many economic traits in some other cultured species [18][19][20]. An increasing number of studies have addressed the transcriptomic profiles of tissues and organs to examine the variety of functional pathways underlying inter-individual difference in feed efficiency [21][22][23][24][25], which focused on target organs such as intestine [21] and liver [22,23], but no similar analysis has been done in neural organs. Moreover, genetic studies on feed efficiency have been scarcely reported in aquaculture species [26,27].
Crucian carp (Carassius auratus) is an economically important aquaculture species worldwide [28]. For this fish, progress on improving feed efficiency has mainly focused on changes in external conditions [29,30]. To date, the only one study about OTL mapping for feed conversion efficiency (FCE, the inverse of FCR) was reported by our laboratory in crucian carp [31]. In a previous study, seven candidate genes associating with FCE and its relevant traits were identified, such as mapk11, cse1l, fam126b, myh14, rgs9bp, cldn10a and cldn10b, involved in several biological functions. So far, few transcriptome analyses have been conducted in crucian carp. Liao, et al. conducted differential gene expression analyses among four tissues [32], and Li, et al. performed RNA-Seq to determine functional differences and DEGs between gibel carps and crucian carp using ten tissues of each species [33]. No transcriptome analysis has been conducted for FCE in this fish. Here, we used high-throughput sequencing of mRNAs from brain tissues of six crucian carps with extreme phenotypes of FCE. We aimed to find novel genes and biological pathways that may be related to feed efficiency of crucian carp, which would provide a better understanding of genetic mechanisms underlying FCE in fish.

Transcriptome Sequencing and Statistics of Unigenes
After quality filtering, the RNA-Seq of six brain samples yielded around 52.59 million high-quality sequence data. The Q30 values of each sample were up to 93.20%, and GC-content of each sample ranged from 45.33% to 49.44% (Table 1). The clean reads obtained from the six transcriptome libraries were assembled to full-length transcripts, and a total of 544,612 unigenes were achieved after elimination of redundant transcripts. Summary data of the assembled transcripts and unigenes are given in Table 2, and the number of unigenes per size that normally decreases continuously is shown in Table 2 and Figure S1.
GO assignment was performed to classify functions of the predicted crucian carp genes. Based on sequence homology, a total of 27,785 unigenes were annotated to 62 terms of GO classification, which were classified into three major functional categories: biological process (21,797 unigenes, 78.45%), cellular component (16,807 unigenes, 60.49%) and molecular function (22,584 unigenes, 81.28%) (Figure 1).
GO assignment was performed to classify functions of the predicted crucian carp genes. Based on sequence homology, a total of 27,785 unigenes were annotated to 62 terms of GO classification, which were classified into three major functional categories: biological process (21,797 unigenes, 78.45%), cellular component (16,807 unigenes, 60.49%) and molecular function (22,584 unigenes, 81.28%) ( Figure 1).

Identification of Differentially Expressed Genes (DEGs)
Overall, 61.15% ± 3.26% of clean reads per sample were mapped back to the assembled transcripts (Table 1); 246 DEGs (Table S1) were identified in the comparison of Low_vs_High FCE groups, with 132 up-regulated and 114 down-regulated genes in the High group, respectively. Among these DEGs, a total of 113 genes were annotated in at least one of the four (Nr, Swissprot, COG and KEGG) databases.

Identification of Differentially Expressed Genes (DEGs)
Overall, 61.15% ± 3.26% of clean reads per sample were mapped back to the assembled transcripts (Table 1); 246 DEGs (Table S1) were identified in the comparison of Low_vs_High FCE groups, with 132 up-regulated and 114 down-regulated genes in the High group, respectively. Among these DEGs, a total of 113 genes were annotated in at least one of the four (Nr, Swissprot, COG and KEGG) databases.

Enrichment for Functional Analysis of DEGs
The results of COG enrichment analysis for the annotated DEGs are shown in Figure 2. Replication, recombination and repair (L, 22-51.16%) annotated the most DEGs in the comparison of Low_vs_High FCE groups, followed by transcription (K, 4-9.3%) and signal transduction mechanisms (T, 4-9.3%).
Pathway analyses were performed using mapped DEGs to other species' orthologues to get more insights into the effects of FCE on metabolisms in the crucian carp. 20 statistically enriched pathways (p < 0.05) were revealed by KEGG pathway analysis, including four in cellular processes, seven in environmental information processing, three in genetic information processing, six in metabolism, and two in organismal systems categories, respectively (Table 3).

Enrichment for Functional Analysis of DEGs
The results of COG enrichment analysis for the annotated DEGs are shown in Figure 2. Replication, recombination and repair (L, 22-51.16%) annotated the most DEGs in the comparison of Low_vs_High FCE groups, followed by transcription (K, 4-9.3%) and signal transduction mechanisms (T, 4-9.3%).
Pathway analyses were performed using mapped DEGs to other species' orthologues to get more insights into the effects of FCE on metabolisms in the crucian carp. 20 statistically enriched pathways (p < 0.05) were revealed by KEGG pathway analysis, including four in cellular processes, seven in environmental information processing, three in genetic information processing, six in metabolism, and two in organismal systems categories, respectively (Table 3).

Critical DEGs Involved in Feed Conversion Efficiency of Crucian Carp
Among all DEGs identified by RNA-Seq analysis, most of the annotated DEGs were classified into molecular function associating with ATP/GTP binding (GO:0005524/GO:0005525), ATP catabolic process (GO:0006200), signal transduction (GO:0007165) and metabolic process (GO:0008152), suggesting that FCE in crucian carp was mainly associated by such kind of pathways as energy metabolism and signal transduction. KEGG analysis also classified DEGs mainly into metabolism and signal transduction, suggesting that these genes may be considered as key candidate genes associating with FCE. Three candidate genes were enriched in metabolism pathways, including Dgkk (Diacylglycerol kinase kappa), Mgst3 (Microsomal glutathione S-transferase 3) and Guk1b (Guanylate kinase 1b). Whereas five genes were classified into signal transduction pathways, such as Vdnccsa1b (Voltage-dependent N-type calcium channel subunit alpha-1B-like) and Tacr2 (Tachykinin receptor 2) in calcium signaling pathway, Tgfα (Transforming growth factor alpha) in ErbB signaling pathway, Vdnccsa1b and Nr4a1 (Nuclear receptor subfamily 4 group A member 1) in MAPK signaling pathway, and Dgkk in phosphatidylinositol signaling system. Attentions to the DEGs with functions involved in growth, such as Endog (Endonuclease G), Crebrtc2 (cAMP-response element binding protein CREB-regulated transcription coactivator 2-like), Myh7 (myosin heavy chain 7), Myh1 (Myosin heavy chain 1), Myh14 (Myosin heavy chain 14) and Igfbp7 (Insulin-like growth factor-binding protein 7), should be also paid because this trait is highly correlated with FCE.

Critical DEGs Involved in Feed Conversion Efficiency of Crucian Carp
Among all DEGs identified by RNA-Seq analysis, most of the annotated DEGs were classified into molecular function associating with ATP/GTP binding (GO:0005524/GO:0005525), ATP catabolic process (GO:0006200), signal transduction (GO:0007165) and metabolic process (GO:0008152), suggesting that FCE in crucian carp was mainly associated by such kind of pathways as energy metabolism and signal transduction. KEGG analysis also classified DEGs mainly into metabolism and signal transduction, suggesting that these genes may be considered as key candidate genes associating with FCE. Three candidate genes were enriched in metabolism pathways, including Dgkk (Diacylglycerol kinase kappa), Mgst3 (Microsomal glutathione S-transferase 3) and Guk1b (Guanylate kinase 1b). Whereas five genes were classified into signal transduction pathways, such as Vdnccsa1b (Voltage-dependent N-type calcium channel subunit alpha-1B-like) and Tacr2 (Tachykinin receptor 2) in calcium signaling pathway, Tgfα (Transforming growth factor alpha) in ErbB signaling pathway, Vdnccsa1b and Nr4a1 (Nuclear receptor subfamily 4 group A member 1) in MAPK signaling pathway, and Dgkk in phosphatidylinositol signaling system. Attentions to the DEGs with functions involved in growth, such as Endog (Endonuclease G), Crebrtc2 (cAMP-response element binding protein CREBregulated transcription coactivator 2-like), Myh7 (myosin heavy chain 7), Myh1 (Myosin heavy chain 1), Myh14 (Myosin heavy chain 14) and Igfbp7 (Insulin-like growth factor-binding protein 7), should be also paid because this trait is highly correlated with FCE.

Discussion
Efficient production of animals plays an important role in livestock [38], poultry [21] and aquaculture [39] industries. This implies the reduction of input costs and the increase of financial returns for producers. In almost all animal production systems, feed cost is the major component (30-70%) of the total costs [40]. Improvement in feed efficiency means reducing the amount of feed resources needed to produce meat and contributing to environmental sustainability [41]. Moreover, it has been reported that feed efficiency is involved in numerous biological processes and functional pathways [10,16], and this trait can be improved by breeding and feeding strategies [24]. Eight QTL associated with FCE, an assessment criterion for feed efficiency, were mapped to four linkage groups, explaining 15.2-20.9% of the phenotypic variations, and seven candidate genes related to FCE and its relevant traits involved in several biological processes were identified in crucian carp in a previous study [34]. In this study, we performed a pioneering and comparative transcriptome analysis based on brain tissues of six crucian carp with extreme FCE phenotypes to help us to understand the molecular mechanisms of FCE. This study is, to the best of our knowledge, the first report of identification of functional determinants involved in FCE at the gene level by transcriptome analysis in fish.
The current RNA-Seq work provided high-quality sequences, and the proportions of mapped reads per sample ranged from 57.73% to 65.35% (Table 1), which ensured the accuracy and reliability of subsequent differential expression analysis. Feed efficiency could be influenced by feed intake to a large content [42], and the brain is a major part of the neural organs that regulates feed intake activity and metabolism in earlier stage of the whole physiological processes in fish. Therefore, fish with significantly different FCE traits may be reflected by RNA profiles in brain. In this study, a total of 544,612 unigenes with a mean size of 644.38 bp (Table 2) were obtained from two groups with extreme FCE performance (namely Low and High groups) of brain tissues in crucian carp, and 246 DEGs (Table S1) that may be involved in FCE were identified. Zhao, et al. identified 300 significantly differentially expressed transcripts in liver of pig [22], while Yi, et al. found 41 promising candidate genes in duodenal of chicken [21]. Compared to these similar reports conducted in other animals, moderate numbers of DEGs were found from brains of crucian carp in this investigation.
According to the results of our previous work conducted using QTL mapping methods [31], identified candidate genes associating with FCE could be classified into three categories: genes encoding GTPase/ATP binding protein, tight junction protein and signal transduction regulator. Moreover, the significantly enriched GO terms of all the predicted crucian carp DEGs in the brain identified in this study were mainly associated with molecular function related to ATP/GTP binding (GO:0005524/GO:0005525), ATP catabolic process (GO:0006200), signal transduction (GO:0007165) and metabolic process (GO:0008152), supporting the hypothesis that FCE of crucian carp are mainly regulated by such pathways as energy metabolism and signal transduction. It has also been reported that energy and vitamin A metabolism pathways in the liver were important for feed efficiency in pigs [22]. Due to relatively a small number of DEGs being annotated in curcian carp, only three candidate genes (Dgkk, Mgst3 and Guk1b) were enriched in KEGG pathways of metabolism, including lipid metabolism, glutathione metabolism, purine metabolism, xenobiotics biodegradation and metabolism that are involved in immune system (Table 3). Dgkk is a member of Dgk (Diacylglycerol kinase) family, which plays an important role in modulating the balance between diacylglycerol and phosphatidic acid [43]. Mgst3 is related to defense mechanisms associated with oxidative stress by utilizing reduced glutathione [44]. For Guk1b, it can affect the metabolism of GTP [45]. A lipid metabolism pathway was the mostly reported pathway that may associate with feed efficiency in other animals [21,24,25,46], suggesting that Dgkk may be considered as an especially critical candidate gene affecting FCE in crucian carp. Moreover, genes related to the immune system were also identified in the transcriptome analyses of feed efficiency in beef [23] and pigs [24].
Another large class of concerned pathways involved in FCE is signal transduction. In this study, calcium signaling pathway (Vdnccsa1b and Tacr2), ErbB signaling pathway (Tgf α), MAPK signaling pathway (Vdnccsa1b and Nr4a1) and phosphatidylinositol signaling system (Dgkk) were found (Table 3) from brain tissues of crucian carp. Among these pathways, ErbB signaling could regulate cell proliferation, differentiation and apoptosis through Akt, MAPK and other signaling pathways. Genes enriched in MAPK signaling pathway related to feed efficiency were generally reported by previous researchers [46,47], such as Fgf (Fibroblast growth factor), Tgfbr (Transforming growth factor beta receptor), Pkc (Protein kinase C) and Hsp72 (72-kd heat shock protein) in beef cattle. Additionally, similar candidate genes associating with FCE have also been found in other animals, which suggested that techniques and strategies employed in these studies were reliable and feasible. Tgfα is an up-regulating gene in High group of crucian carp involved in ErbB signaling pathway in this study, while Tgfbr is also associated with feed efficiency in chicken [48] and common carp [26]. Igf1 (Insulin-like growth factor 1) is an important regulator of muscle growth and energy metabolism [49,50] that was identified as a candidate gene involved in feed efficiency in chicken [51] and common carp [26]. Igfbp3 (Insulin-like growth factor-binding protein 3) and Igfbp7 are Igf binding proteins in circulation and played important roles in modulating Igf bioavailability and half-life, which were respectively found in beef [23] and crucian carp. Both Myh7 and Myh1 were digged out in this report that may be involved in regulating muscle growth, and Myh7 has also been found to be differentially expressed between low-and high-FCR groups in pigs [22]. Moreover, Myh14, another gene from Myh superfamily that plays an important role in fish muscle formation [52], has been identified in our previous work conducted using QTL mapping method in crucian carp [31]. The consistent result obtained in both studies suggesting that Myh14 may be also considered as an especially critical candidate gene affecting FCE in crucian carp.
The reliability of RNA-Seq analysis was confirmed by the concordance between the computational and experimental results by employing nine randomly selected DEGs for qRT-PCR assays in this study, which was similar to some previous studies in animals [53][54][55]. Since there is no reference genome for Carassius auratus, most unigenes were not well annotated in all annotation databases. At present, a well-assembled genome of crucian carp is needed, and it may help to find more key candidate genes potentially related to feed efficiency. Our results would provide useful information for a physiological basis to develop improved feed formulas and/or feeding conversions. Furthermore, these novel findings would be useful for such future studies as association analysis and functional verification between candidate genes and FCE traits, thereafter providing proper markers for marker-assisted selection for potential improvement of feed efficiency in crucian carp.

Sample Collection and RNA Preparation
A reference family consisting of 120 diploid crucian carp fingerlings, generated by parents from Zhangdu Lake, Yangtze River (Wuhan, China) using artificial crossing in April, 2015, was used in this study. At 82 days post hatch (dph), these reference fish (mean body weight of 0.87 ± 0.39 g) were individually reared in the same environment with pallet feed in a series of re-circulating aquarium tanks in order to achieve accurate feed consumption. Such detailed conditions of the aquarium tanks as water temperature (27-28 • C), dissolved oxygen (7-8 mg/L) and water flow rate (1 ms −1 ), were regularly maintained throughout the feeding test for two months. During the experimental period, all test fish were fed three times (10:00 am, 15:00 pm and 20:00 pm) a day by the same fish-feeder, and the fish-feeder stopped feeding when the fish satiated each meal. Faeces generated by experimental fish in each tank were siphoned out once a day and a complete water change was made once a week. The individual body weight (BW) was recorded at the start (initial BW, BW I ) and the end (final BW, BW F ) of the experiment to calculate BW gain. The feed conversion efficiency (FCE) was estimated with the model as follows: Here FCE = feed conversion efficiency, BW F = final body weight, BW I = initial body weight, FI = total feed intake of each individual, which was recorded as the difference between the final and the beginning weight of diet used during the test. The detailed information of this feeding test was described in previous study [31].
Fish with extreme FCE phenotypes were used for RNA extraction. At the end of the whole experimental period, we selected six fish consisting of two groups (three biological replicates per group) to represent two divergent FCE performances (Low group: L1-L3 and High group: H1-H3). Table 4 detailed the measurements of FCE and its relevant traits in two groups. For RNA isolation, brain tissues were sampled immediately from sacrificed fish, frozen in liquid nitrogen and then stored at −80 • C until further processing. All experimental procedures involved in fish in this study were based on institutional regulations and guideline for experimental animals of the Hubei Provincial Committee for Animal Welfare (Permit Number: 20130522-02).

RNA Extraction, Library Preparation and Transcriptome Sequencing
Total RNA was isolated from the frozen brain samples using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The RNA degradation and contamination was monitored on 1% agarose gels. RNA purity and concentration was measured by NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA) and Qubit ® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). The RNA integrity number was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), which ranged from 8.2 to 8.9, suggesting that the samples were well preserved to meet the cDNA library construction requirements.
Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) according to manufacturer's recommendations. Oligo-dT beads (Qiagen, Dusseldorf, Germany) were used to separate poly (A) mRNA from the total RNA of each sample, and the fragmentation buffer was added to split all mRNA into short fragments. Then the first-strand cDNA was synthesized using the random hexamer-primed reverse transcription, and the second-strand cDNA was generated using RNase and DNA polymerase I. The cDNA fragments were washed by EB buffer for end reparation poly (A) after purification, and index codes were added to attribute sequences to each sample. Fragments with suitable size for sequencing were isolated from the agarose gels, and PCR amplification was used to enrich these fragments to construct the final cDNA libraries. Finally, the cDNA libraries were sequenced on an Illumina Hiseq 2000 platform and paired-end reads were generated.

Transcriptome Assembly and Functional Annotation
For ensuring high-quality data, raw data (raw reads) of fastq format were firstly processed through in-house Perl scripts, which eliminated all those reads with sequencing adapter and nucleotides in reads with quality value less than 20 in both end. In this step, good quality sequences of clean data (clean reads) were obtained by abandoning reads containing adapter, reads containing ploy-N and low-quality reads from raw data. At the same time, Q30, GC-content and sequence duplication level of the clean data were calculated. The clean data of this article are publicly available in the NCBI Sequence Reads Archive (SRA) with accession number PRJNA433432.
All clean reads of the six libraries were jointly assembled into contigs employed by Trinity software [56]. Since there is no reference genome of Carassius auratus, a k-mer value cutoff of 25 was used after removing redundant nucleotide sequences by Tgicl (v2.1, http://sourceforge. net/projects/tgicl/files/tgicl%20v2.1/). Then, unigenes were generated by connecting the contigs (longer than 200 bases) to obtain sequences that could not be extended on either end, and maximum length non-redundant unigenes were acquired by further splicing and assembling using TGICL clustering software (J. Craig Venter Institute, Rockville, MD, USA). Finally, unigenes were aligned against the Nr (NCBI non-redundant protein sequences), Swissprot (A manually annotated and reviewed protein sequence database), COG (Clusters of Orthologous Groups of proteins), and KEGG (Kyoto Encyclopedia of Genes and Genomes) of protein databases using BlastX with an E-value <10 −5 . GO (Gene Ontology) annotation of these unigenes was performed using Blast2GO (https://www.blast2go.com/) based on the results of the NCBI Nr database annotation. Blastn was used to align these unigenes with the Nr database to search for proteins with the highest sequence similarity to the given unigenes and annotate their protein functions at the same time.

Analysis of Differentially Expressed Genes (DEGs)
Gene expression levels were estimated by RSEM (RNA-Seq by Expectation Maximization) [57] software package for each sample. The mapped reads were normalized according to fragment per kilobase of exon model per million mapped reads (FPKM) for each unigene between the two groups (Low_vs_High). Differentially expressed genes (DEGs) between the two groups were identified by the DEGseq package (samples with three biological replicates) applying the MA-plot-based method with Random Sampling model (MARS) method. In this study, DEGs with significant expression abundance between the two groups were selected using the following filter criteria: p-value < 0.01 and the absolute value of log2 Ratio ≥ 1, meaning each DEG between two groups should be at least two-fold. In order to determine the potential functions and metabolic pathways of these DEGs, COG, GO and KEGG enrichments were further analyzed. COG annotation of the DEGs was performed using Blastall software. GO enrichment analysis (p-value ≤ 0.05) of the DEGs was implemented by the topGO R packages based Kolmogorov-Smirnov test. Based on the hyper-geometric distribution model, we used KOBAS software [58] to test the statistical enrichment of DEGs in KEGG pathways, and the enrichment p-values were adjusted using the Benjamin and Hochberg method.

Validation of RNA-Seq Results by qRT-PCR
To confirm our differential expression results of RNA-Seq, we conducted quantitative reverse transcription PCR (qRT-PCR) assays for nine randomly selected DEGs in the same RNA samples used for RNA-Seq. First-strand cDNA was reverse-transcribed from total RNA using Reverse Transcriptase M-MLV (TaKaRa, Tokyo, Japan) with oligo-dT primer following the manufacturer's instructions. The first-strand cDNA from each sample was diluted by 1:5 with sterile distilled water and used as template. Primers (Table S2) for qRT-PCR analyses of the nine genes were designed using Primer 5 Software, and qRT-PCR was performed on a StepOneTM Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The qRT-PCR reaction solution consisted of 6.5 µL Power SYBR Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA), 0.2 µM of each forward and reverse primer, 1.2 µL diluted cDNA and 4.5 µL sterile distilled water. PCR program was 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 30 s, 72 • C for 45 s. Three parallel experiments were conducted in each run for each sample. The relative expression levels were normalized towards the internal control gene of β-actin. Optimized comparative Ct value (2 −∆∆CT ) method [59] was used here to estimate gene expression levels.

Conclusions
We have reported a cerebric transcriptome of crucian carp using six fish with extremely low and high FCE. A total of 544,612 unigenes with a mean size of 644.38 bp were obtained from two FCE groups, and 246 DEGs that may be involved in FCE were identified. Based on GO and KEGG annotations, 13 DEGs related to metabolism, signal transduction and growth were identified to be key candidate genes associated with FCE traits. Our results provide valuable information for elucidating molecular mechanisms of feed conversion efficiency in fish, and these novel findings would be useful for such future studies as association analysis and functional genomics verification of candidate genes related to FCE in crucian carp.