Comparative Transcriptomic Study of Muscle Provides New Insights into the Growth Superiority of a Novel Grouper Hybrid

Grouper (Epinephelus spp.) is a group of fish species with great economic importance in Asian countries. A novel hybrid grouper, generated by us and called the Hulong grouper (Hyb), has better growth performance than its parents, E. fuscoguttatus (Efu, ♀) and E. lanceolatus (Ela, ♂). We previously reported that the GH/IGF (growth hormone/insulin-like growth factor) system in the brain and liver contributed to the superior growth of the Hyb. In this study, using transcriptome sequencing (RNA-seq) and quantitative real-time PCR (qRT-PCR), we analyzed RNA expression levels of comprehensive genes in the muscle of the hybrid and its parents. Our data showed that genes involved in glycolysis and calcium signaling in addition to troponins are up-regulated in the Hyb. The results suggested that the activity of the upstream GH/IGF system in the brain and liver, along with the up-regulated glycolytic genes as well as ryanodine receptors (RyRs) and troponins related to the calcium signaling pathway in muscle, led to enhanced growth in the hybrid grouper. Muscle contraction inducing growth could be the major contributor to the growth superiority in our novel hybrid grouper, which may be a common mechanism for hybrid superiority in fishes.


Introduction
Grouper (Epinephelus spp.) is a group of important economic marine fish species that are widely cultured in China and Southeast Asian countries [1,2]. However, some threats, such as germplasm resource degradation and insufficiency of high-quality gametes in hatcheries, restrict the development of grouper aquaculture [3]. One critical way to solve these problems is hybridization technology, which has been widely applied to grouper artificial breeding over the past few decades [4][5][6][7][8][9]. A novel hybrid grouper (Hyb), generated by us and called the Hulong grouper, exhibits a faster growth rate and stronger disease resistance than its parents a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 regulated in the hybrid grouper compared to its parents. In addition, the calcium signaling pathway was also highly activated and the expressions of troponins, including TnC, TnT and TnI were much higher. These findings enhanced our understanding of pathways regulating muscle development and heterosis in fish.

Materials and Methods
Fish and sample preparation E. fuscoguttatus (♀), E. lanceolatus (♂) and their hybrid F1 offspring (the Hulong grouper) were cultivated under the same breeding conditions (i.e., maintained in a laboratory recirculating seawater system at 25-30˚C) in the Daya Bay Seawater Fish Farm in Huizhou, Guangdong Province, China [21]. Three individuals from each species were randomly selected at the age of 18 months. Fresh muscle tissues were dissected from each fish after euthanasia by immersion in MS-222 buffered solution (3 g/L) on ice, and the tissues were immediately deposited in separate sterile tubes soaked in liquid nitrogen. The research protocol and the procedure for handling experimental animals were reviewed and approved by the Institutional Review Board on Bioethics and Biosafety of BGI.

RNA extraction and high-throughput sequencing
Total RNA from nine muscle samples was separately isolated with Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified using the RNeasy Animal Mini Kit (Qiagen, Valencia, CA). RNA quality was then assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, Calif). An equal amount of total RNA with a 28S/18S ratio ! 1.0 and RNA Integrity Number (RIN) ! 7.0 from three individuals in each group was pooled together to construct cDNA libraries (Illumina, San Diego, CA). These libraries were subsequently sequenced through an Illumina HiSeq2000 platform at BGI (Shenzhen, China). Approximately 5 Gb of raw reads were generated for each of the three samples.

Raw reads processing, RPKM calculation and GO enrichment
Adaptor sequences and raw reads with a more than 10% N bases and quality score below 20 were filtered out using SOAPnuke [53]. The remaining clean reads were then aligned to the available genome of E. coioides (from our unpublished data) with SOAPaligner 2 [54]. The Cuffdiff package of Cufflink software (version 2.1.1) [55], with the core parameters (-FDR 0.05-geometric-norm TRUE-compatible-hits-norm TRUE) to reduce certain types of bias caused by differential amounts of RNA reads, was used to calculate the RPKM (Reads Per Kilobase Transcriptome per Million mapped reads) values of the mapped genes. The differential expressed genes (DEGs) were identified using the absolute value of log2 (ratio) !1 as the threshold. Gene ontology (GO) annotation of genes was conducted with Blast2GO [56]. The enrichment analysis was carried out by GO:Termfinder using the hypergeometric test [57,58] and P-values were corrected using the Bonferroni method [59].

Validation of the RNA-seq analysis by qRT-PCR
A qRT-PCR process is typically used to confirm data obtained from high-throughput sequencing [60,61]. In this study, 16 identified genes were used to verify the analysis results of transcriptional sequencing. Total RNA from the muscle tissue of the hybrid grouper and its parents were extracted and purified with the RNeasy Animal Mini Kit (Qiagen, Valencia, CA), and first-strand cDNA was obtained using a RevertAid First Strand cDNA Synthesis Kit (Fermentas, Vilnius, Lithuania). The qRT-PCR reactions based on SYBR (SYBR Green I, Osaka, Japan) were performed with a LightCycler 480 system (Roche, Basel, Switzerland) using the reaction procedure as follows: 95˚C for 1 min followed by 40 cycles of 95˚C for 10 s and 60˚C for 30 s. All samples were examined in triplicate (i.e., three separate biological replicates) with β-actin as the internal control, and the 2 -ΔΔCt method [62] was used to calculate relative expression amounts. A standard curve was generated by running serial 10-fold dilutions (6 × 10 2 to 6 × 10 8 ) of a recombinant plasmid (pLB-β-actin) obtained as the template. Copy numbers of βactin were detected to determine its stable and uniform expression in all samples. All gene-specific primers for target genes as well as the β-actin are summarized in the S1 Table. Results

High-throughput sequencing and mapping to the reference genome
To obtain the muscle transcriptome data of the hybrid F1 and its parent non-hybrids (1.5 years was used as the approximate age for each fish to remove age difference) and to investigate the gene expression variations among the 3 grouper species, three individuals from each species were caught and a total of 9 RNA samples were isolated from muscle tissues. All RNA samples were of high-quality with a 28S/18S ratio ! 1.0 and RIN ! 7.0 after an RNA quality assessment. Three cDNA libraries were then constructed separately for the Hyb, Efu and Ela. Approximately 5 Gb of raw data for each sample were produced, and a total of 169 million paired-end clean reads were finally obtained after data filtration ( Table 1). All transcriptome reads generated in this study have been submitted to the Sequence Read Archive (SRA) database of NCBI (Hyb: SRX1631685; Efu: SRX1626373; Ela: SRX1631646).
The clean reads were subsequently aligned to the E. coioides genome (from our unpublished data) with SOAP aligner 2.0 [54]. We assembled the reference based on more than 100× coverage paired-end reads from libraries of 200-bp, 500-bp, 800-bp, 2-kb, 5-kb, 10-kb and 20-kb insert sizes to achieve a good genome assembly. Approximately 62.43% (Hyb), 60.93% (Ela) and 60.79% (Efu) of clean reads were mapped to the reference genome ( 5 base mismatches), respectively, in which 54.06% (Hyb), 51.05% (Ela) and 55.71% (Efu) of the reads were mapped to the gene regions ( 5 base mismatches). Moreover, 53.27% (Hyb), 48.59% (Ela) and 52.09% (Efu) of those uniquely mapped reads were further selected for gene quantification analysis (Table 1). Interestingly, all these mapping rates are higher than the previously published data on the liver and brain [21], which suggested that genes expressed in muscle may be more conserved than those in the other two tissues among the studied groupers.

RPKM calculation and GO enrichment of DEGs
Gene expression levels were quantified by RPKM values. RPKMs of each gene in the muscle of the hybrid grouper were compared to its parents (S2 Table). The differentially expressed genes (DEGs) between the hybrid grouper and its parents were identified through filtering based on the criteria of false discovery rate (FDR) 0.05 and absolute log2 (ratio) !1. Our data showed that compared to Efu, there were 553 up-regulated and 1,250 down-regulated genes in the Hyb (S3 Table). Compared to Ela, the corresponding numbers of DEGs were 389 and 3,287, respectively (S4 Table). In total, there were 5,479 genes that were differentially expressed in the Hyb compared to its parents, and among those genes, 574 overlap DEGs were found in 'Hyb vs. Efu' and 'Hyb vs. Ela'. Therefore, a total of 4,905 unique DEGs were identified and they were further selected for Gene Ontology (GO) enrichment analysis. GO analysis of our dataset by Blast2GO [56] showed that all the above-mentioned 4,905 unique DEGs identified in the Hyb were enriched into 55 GO terms, in which cell and cell part were the top two enriched terms with 1,357 and 1,356 annotated genes, respectively (Fig 1). Among these 4,905 genes, 1,715 genes were grouped into biological process, 1,698 genes into cellular component and 1,771 genes into molecular function (Fig 1). In the liver, we previously found that metabolic process and catalytic activity were the most enriched processes out of the molecular function and biological process categories [21]. In the muscle, a large number of genes were also grouped into metabolic process (868 genes) and catalytic activity (1,050 genes), which made them the second and third most-enriched terms in each categories. The results suggested that extensive metabolic and catalytic activities in the muscles of the hybrid may also contribute to enhanced growth.

Glycolytic genes are up-regulated in the hybrid grouper
We formerly showed that in the liver and brain the upstream GH/IGF axis and its downstream signaling pathways, such as glycogen synthesis, may be related to growth superiority in the hybrid grouper [21]. In this study, we further investigated whether in muscle, as one of the major tissues affected by the GH/IGF pathway, the genes in related downstream pathways (such as glycolysis) expressed differently between the hybrid and its parents.
Glycolysis consists of 10 major steps with each reaction catalyzed by a specific enzyme. RPKM values of these 10 enzymes (Pgm, Gpi, Pfk, Ald, Tpi, Gapdh, Pgk, Pgam, Eno, Pyk) in the muscle of Hyb were compared to RPKM values of its parents Efu and Ela (Table 2). In the Hyb, we observed that all the 10 candidate genes had higher RPKMs than Ela with an excess number ranging from 280.91 (Pfk) to as high as 15,450.35 (Ald). Compared to Efu, up to 8 (Pgm, Gpi, Pfk, Ald, Gapdh, Pgk, Pgam, Eno) out of the 10 glycolytic genes were up-regulated in the Hyb, with an excess number ranging from 41.68 (Pgk) to as high as 7,002.71 (Gapdh). Consistent with this outcome, the RPKM of Ldh, which codes for the protein catalyzing the conversion of lactate to pyruvic acid, was much higher (6,245.80) in the Hyb than in the Efu (4,967.76) and Ela (4,986.57), which provides further evidence that the glycolysis pathway was activated in the Hyb (Fig 2). These genes may lead to active anabolism and catabolism in muscle of the Hyb and subsequently stimulate the development of the muscle in the hybrid grouper [34,63].

Ca 2+ signaling is activated and Tns were up-regulated in the Hyb
We further investigated, in addition to glycolysis, if genes in other pathways involved in muscle growth demonstrated different expression levels between the Hyb and its parents as well. Troponin (Tn) is a complex of three regulatory proteins (troponin C, troponin T and troponin I) that is integral to muscle contraction and development [46,47]. Troponin C (TnC), containing four Ca 2+ -binding EF hands [42,43], can directly bind to and sense Ca 2+ [37,44,45], which was released into the cytosol by RyRs cellular calcium channels [38,39]. We thus tested whether troponins and genes involved in the calcium signaling pathway were also up-regulated in the hybrid grouper.
By comparing the RPKMs of Ca 2+ -regulating protein-coding genes in the Hyb and its parents ( Table 3, Fig 3), RyRs (RyR1 and RyR3) were found to be up-regulated in the muscle of Hyb, which suggested there were active changes in the Ca 2+ concentration of the cytosol in the cytoplasm. The downstream Ca 2+ -binding protein TnC was also up-regulated in the Hyb along with the other two troponins, TnI (TnI1, TnI2) and TnT (TnT1, TnT2) ( Table 3, Fig 3). The activation of Ca 2+ signaling as well as the abundance of the Tn complex may lead to severe muscle movement which facilitated the increase in the size of skeletal muscle [35].

Validation of the RNA-seq analysis by qRT-PCR
To further verify the expression levels of the above-mentioned genes in glycolysis and Ca 2+ signaling pathways, qRT-PCR validation experiments were conducted using gene-specific primers (S1 Table), which were designed based on the mapped genome sequences. The efficiency of the primers was determined using PCR. A single band (ranging from approximately 150 bp to 300 bp in length) was observed in each lane, and the PCR products were sequenced and analyzed. The result showed that they were the target gene fragments. Standard curve analysis of β-actin showed R2 value of 0.99 and slope of -3.48, from which E (amplification efficiency)  was calculated to be 94%, which is acceptable for qRT-PCR analysis. The calculation of copy numbers showed that β-actin exhibited stable and uniform expression in the three species, suggesting that it could be used as the housekeeping gene in this study. The data showed that, consistent with the RNA-seq analysis, all glycolytic genes had higher expression in the Hyb than in Ela. Compared to Efu, Pgm, Pfk, Gapdh, Pgk and Eno were upregulated in the hybrid grouper (Fig 4A). Regarding the troponins and genes involved in Ca 2+ signaling pathway, all those genes had higher expressions in the Hyb compared to the Ela, whereas all troponins were up-regulated in the muscle of the Hyb compared to Efu (Fig 4B). These results confirmed that the expression differences of most studied transcripts among the 3 groupers revealed by both methods were generally consistent.

Discussion
A novel grouper hybrid called the Hulong grouper (Hyb) exhibits better growth performance compared to its parents, Ela and Efu [10,11]. In a previous study, we had inferred that the GH/IGF system and its downstream signaling pathways, such as glycogen synthesis, may contribute to growth superiority [21]. However, the underlying mechanism for muscle growth superiority in the Hyb is still unknown. Here, transcriptome analysis of muscle in Hulong grouper and its parents had been carried out, and we identified candidate genes regulating downstream glycolysis as well as muscle contraction involved Ca 2+ signaling pathway and troponin synthesis. The qRT-PCR confirmed the differential expression in the selected genes. All of these differences, to a great degree, may account for the enhanced muscle growth of the hybrid.

Glycolysis involved in muscle development
Glycolysis was reported to be involved in muscle development. Mutation of glycolytic genes, such as Pfkm, Pgam and Pyk, leads to different disease symptoms in muscle [34,64]. Multiple muscle diseases have been known to show a reduction in glycolytic enzyme expression and an accumulation of glycogen [29][30][31][32][33]. In addition, reduced activity of glycolytic enzymes were observed in experimental animals with muscle-specific inactivation of proteins (such as mTOR) regulating cell growth [65], while significantly up-regulated expressions of glycolytic genes were detected in skeletal muscle-specific, conditional transgenic mice that expressed extra exogenic proteins (such as Akt1) stimulating muscle fiber growth [66]. In this study, almost all glycolytic genes were up-regulated in the muscle of the Hyb compared to its parents, which indicated that glycolysis may play an important role in the superior growth of the hybrid grouper.
Many metabolites in glycolysis are intermediates for other anabolic pathways as well, and as a consequence, the activation of the glycolysis pathway also led to the activation of multiple metabolic pathways, such as gluconeogenesis [67,68], lipid metabolism [69], the citric acid cycle (CAC, also known as tricarboxylic acid or TCA cycle, a post-glycolytic process leading to amino acid synthesis, nucleotide synthesis and tetrapyrrole synthesis) [70,71], among other processes. Therefore, the up-regulated expression of glycolytic enzymes may synchronously influence growth superiority in the Hyb not just via glycolysis, but in a much wider scope. Glycolysis is a downstream pathway of the GH/IGF axis [72,73]. Our previous study suggested that variation of GH/IGF system in the brain and liver may be the important cause of enhanced growth in the Hyb [21]. In our current study, we confirmed that in muscle, a target of the GH/IGF system, the downstream glycolysis was also activated in the Hyb with up-regulation of almost all glycolytic genes, which provided evidence for a possible critical role in the upstream GH/IGF pathway in muscle development by stimulating downstream glycolysis in the muscle as the major effector organ (Fig 5).

Ca 2+ signaling pathway and calcium-binding protein TnC
Calcium acts as a signal messenger combining excitation events with downstream effects [74]. Muscle contraction results from the binding of Ca 2+ to a calcium-binding protein TnC in skeletal muscle [75,76]. TnC, in addition to TnI and TnT, forms the troponin complex that binds were picked out as representatives of differences between the Hyb and its parents. β-actin was used as the internal control, and each value represents an average of three separate biological replicates. The error bars represent the standard deviation. Asterisk (*) was marked if there was significant difference between the two samples, as assessed by a t-test (P < 0.05). doi:10.1371/journal.pone.0168802.g004 to actin and tropomyosin [75]. The release of Ca 2+ into the cytosol from SR via the RyR receptor channels raised the Ca 2+ concentration and its binding to TnC changes the conformation of TnC, which ultimately resulted in the strong interaction between actin and myosin to generate force and muscle shortening [77]. In the muscle of the Hyb, Ca 2+ signaling was activated and troponins were up-regulated, which suggested there was an active state of muscle contraction that led to fiber growth [35]. These data were consistent with the active swimming of the Hyb.

An overall view on the the superior growth of the Hyb
The heterosis of the hybrid Hyb grouper is represented in many aspects, such as faster growth rate and stronger disease resistance [10,11]. Although the biological and physiological features of the Hyb have been well-studied, the genetic mechanisms of its heterosis, especially enhanced growth, had yet to be investigated. We formerly obtained the transcriptome sequences of the tissues in the neuroendocrine system (brain and liver) from the Hyb and its parents. Differential gene expression was observed in the GH/IGF axis and downstream signaling pathways, such as protein synthesis and glycogen synthesis [21]. In this study, we sequenced the transcriptomes of the major effector organ (muscle) from these three groupers, and the downstream glycolysis pathway was shown to be highly activated in the Hyb. In addition, the Ca 2+ signaling pathway was activated via up-regulated RyRs (RyR1 and RyR3) and troponins (TnC, TnI and TnT), which resulted in elevation of muscle contraction and subsequently promoted muscle growth and body development (Fig 6).

Conclusions
The novel hybrid grouper, called the Hulong grouper (Hyb), exhibits better growth performance and stronger disease resistance than its parents Efu (♀) and Ela (♂). Using the RNA-seq technique, we obtained the muscle transcripts from these 3 species. The transcriptomic data showed that almost all glycolytic enzymes were up-regulated in the Hyb. Troponins as well as genes involved in Ca 2+ signaling were also up-regulated. The RNA expression levels of candidate genes were further validated by qRT-PCR. These findings revealed that glycolysis in addition to calcium signaling and up-regulated troponins may be the most important contributors to the growth superiority of the hybrid grouper.
Supporting Information S1