Transcriptome Profiling Revealed Basis for Growth Heterosis in Hybrid Tilapia ( Oreochromis niloticus ♀ × O. aureus ♂ )

: Hybrid tilapia were produced from hybridization of Nile tilapia ( Oreochromis niloticus ) and blue tilapia ( O. aureus ). Comparative transcriptome analysis was carried out on the liver of hy ‐ brid tilapia and their parents by RNA sequencing. A total of 2319 differentially expressed genes (DEGs) were identified. Trend co ‐ expression analysis showed that non ‐ additive gene expression accounted for 67.1% of all DEGs. Gene Ontology and KEGG enrichment analyses classified the re ‐ spective DEGs. Gene functional enrichment analysis indicated that most up ‐ regulated genes, such as FASN, ACSL 1 , ACSL 3 , ACSL 6 , ACACA, ELOVL 6 , G 6 PD, ENO 1 , GATM , and ME 3 , were involved in metabolism, including fatty acid biosynthesis, unsaturated fatty acid biosynthesis, glycolysis, pentose phosphate pathway, amino acid metabolism, pyruvate metabolism, and the tricarboxylic acid cycle. The expression levels of a gene related to ribosomal biosynthesis in eukaryotes, GSH ‐ Px , and those associated with heat shock proteins (HSPs), such as HSPA5 and HSP70 , were significantly down ‐ regulated compared with the parent tilapia lineages. The results revealed that the metabolic pathway in hybrid tilapia was up ‐ regulated, with significantly improved fatty acid metabolism and carbon metabolism, whereas ribosome biosynthesis in eukaryotes and basal defense response were significantly down ‐ regulated. These findings provide new insights into our understanding of growth heterosis in hybrid tilapia.

China is the largest producer of tilapia with a production about 1.65 million tons in 2020 [22]. Crossbreeding is a common strategy used for genetic improvement within tilapia species [23,24]. Hybrid tilapia, produced from the hybridization of maternal O. niloticus and paternal O. aureus, are popular due to the high male percentage, which contributed significantly to the total tilapia production in China [24]. Hybrid tilapia shows a faster growth rate when compared to its parents [14]. However, the molecular bases are not well understood. Some authors have highlighted anthropogenic effects [25,26]. Theoretically, no new genes are produced in the F1 hybrid, so heterosis is likely caused by differences resulting from qualitative or quantitative modifications of gene expression [27]. Gene expression patterns can be divided into either additive or non-additive expression on the basis of differential expression in the hybrids compared with the parents. Nonadditive genes, which are genes in the hybrids that show significantly different expression from the mid-parent value, have been suggested to be associated with heterosis [28][29][30]. Recent studies of hybrid yellow catfish Huangyou-1 (Pelteobagrus fulvidraco ♀ × P.vachelli ♂) [8] and hybrid abalone (Haliotis diversicolor Taiwan genotype ♀ × H. diversicolor Japan genotype ♂) [31] have also focused on the non-additively expressed genes to characterize gene regulation and its underlying mechanisms.
In order to reveal the mechanisms underlying the heterosis in hybrid tilapia (O. niloticus ♀ × O. aureus ♂), comparative transcriptome analysis was carried out between hybrid tilapia and its parental lines. It is hypothesized that the non-additive modes of gene expression are vital to realization of heterosis generation in hybrid tilapia.

Sample Collection and Preparation
Experimental fish used in the study were cultivated at the Tilapia Genetic Breeding Center of the Ministry of Agriculture and Rural Affairs, Wuxi, China. Two tilapia varieties, "Egypt strain" O. niloticus (N for short) and "Xia'ao strain" O. aureus (A for short), were both from inbred lines produced through several generations of sibling mating; the sterile hybrid tilapia (N ♀ × A ♂, NA for short) were the F1 offspring of these two lines. All individual fish were cultured in three open 0.067-ha pools under the same conditions. Then one-hundred healthy fish of each line were randomly selected and kept in a recirculating system (360-L tanks with circulating aerated water) to acclimate for 14 days. Water quality was monitored daily and maintained as follows: temperature, 29.0 ± 1 °C; pH, 7.6-7.8; dissolved oxygen (DO), >5.0 mg/L. The fish were fed a commercial feed twice a day (8:00 and 16:00). No tilapia died during the acclimation.
At the end of the trial, after 24 h of starvation, 30 males of each line were randomly selected for measure. Hybrid tilapia had a body mass of 63.65 ± 4.84 g. The body mass of O. niloticus was 59.65 ± 4.59 g, whereas that of O. aureus was 60.71 ± 5.63 g, which were lighter than that of hybrid tilapia (p < 0.05). Three males per line with a normal appearance and good vitality were deeply anesthetized with 100 mg/L MS-222 (Sigma-Aldrich, St. Louis, MO, USA) before dissection. There were no skin lesions before tissues were collected. The liver tissues excised from fish were snap-frozen in liquid nitrogen and then stored at −80 °C for subsequent RNA isolation. The protocols used meet the Zhejiang Laboratory Animal Management guidelines (000014349/2017-768264) established by the government of Zhejiang province on the use and care of animals.

RNA Extraction, Library Construction, and Transcriptome Sequencing
Total RNA was extracted using TRIzol™ Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Three independent biological replicates were performed for each tilapia line. RNA concentration was measured using a NanoDrop™ 2000 (Thermo Scientific, Waltham, MA, USA). RNA integrity was assessed using the RNA 6000 Nano Assay Kit of the Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA).
For RNA library preparation, the RNA of each individual from each tilapia line was diluted to the same concentration with RNase-free water. The three RNA samples from each line were equally mixed together (1 μg RNA from each sample) to obtain three groups. In order to select cDNA fragments of preferentially 240 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Sequencing libraries were generated using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumina, San Diego, CA, USA). After cluster generation, the library preparations were sequenced on an Illumina HiSeq X Ten platform by Biomarker Technologies (Beijing, China). Paired-end reads were generated, each about 150 bp long.

Data Filtering, Comparative Analysis, and Gene Functional Annotation
Raw reads in FASTQ format were first processed through in-house Perl scripts. In this step, clean reads were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads from the raw data. Simultaneously, Q30, GC content, and sequence duplication level of the clean data were calculated. The clean data were submitted to the National Center for Biotechnology Information (NCBI), deposited in the Sequence Read Archive (SRA) under BioProject PRJNA792260(SUB10844155), with the accession numbers SRR17399951, SRR17399952 and SRR17399953.
These clean reads were then mapped to the reference genome sequence (https://www.ncbi.nlm.nih.gov/genome/197?genome_assembly_id=28435, accessed on 18 March 2015). Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. The TopHat2 tool was used to map the reads to the reference genome [32].
Gene function was annotated based on the following databases: NCBI non-redundant protein sequences, NCBI non-redundant nucleotide sequences, Pfam (Protein family), Clusters of Orthologous Groups of proteins, Swiss-Prot, KEGG Ortholog, and Gene Ontology (GO).

Quantification of Gene Expression Levels and Differential Expression Analysis
Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). Differential expression analysis between hybrid tilapia and their parents was performed using the EBSeq R package [33]. The resulting p values were adjusted using the Benjamini and Hochberg approach to control the false discovery rate. Genes showing a false discovery rate < 0.01 and a |log2(fold-change)| ≥ 1 were considered significant differentially expressed genes (DEGs). The difference in gene expression level between two tilapia lines was displayed by volcano plot. The DEGs were subjected to trend co-expression analysis according to K-means clustering and were classified into five patterns ( Figure 1): over-dominant, high-parent dominant, low-parent dominant, under-dominant, and additive. These analyses were performed using BMKCloud (www.biocloud.net, 7 March 2020).

GO and KEGG Enrichment Analysis of DEGs
GO enrichment analysis of DEGs was implemented in the goseq R package using the Wallenius non-central hyper-geometric distribution [34]. KOBAS 2.0 (Peking University, Beijing, China) was used to test the statistical enrichment of DEGs in KEGG pathways [35]. GO terms and KEGG pathways showing a corrected p ≤ 0.05 were considered significantly enriched. Hierarchical cluster analyses of DEGs were constructed. The R program was used to depict the heat maps for gene clustering.

Validation of RNA-Seq Data by Quantitative Real-Time PCR (qRT-PCR)
Sixteen growth-related genes were selected for RNA-Seq data validation using qRT-PCR. After transcriptome sequencing, the remaining nine RNA samples were returned by mail on dry ice; then, 1 μg of total RNA was converted to cDNA by a PrimeScript™ RT reverse transcription kit (TaKaRa, Dalian, China) according to the manufacturer's instructions. Primers were designed using Primer Premier 5 and were synthesized at Biotech Bioengineering Co. (Shanghai, China) ( Table 1)

Statistical Analysis
All qRT-PCR data were expressed as the mean ± standard error. Statistical significance of differences was determined using one-way ANOVA in SPSS 25.0 (IBM, Chicago, IL, USA). p < 0.05 was considered to indicate a statistically significant difference. Tukey's test was applied to detect significant differences (p < 0.05). Graphs were made with GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA).

Sequencing Data Statistics
After sequence quality control, a total of 19.54 Gb clean data were obtained, and the Q30 base percentage was above 90.20% (Supplementary Table S1). RNA-Seq obtained over 6.22 Gb of clean reads from each sample, and the average GC content was approximately 50%. The alignment efficiency ranged from 65.77% to 73.85% (Supplementary Table S2). A total of 1687 novel genes were discovered, of which 1546 were annotated.  Figure S1).

Identification of Differentially Expressed Genes
Of the 1440 up-regulated DEGs between the pure species (N and A) and hybrid (NA), 220 DEGs were also among the DEGs between N and A. In addition, 163 DEGs were specifically expressed in "N vs. NA" and 420 DEGs in "A vs. NA"; 4 DEGs were co-expressed in "N vs. NA", "A vs. NA", and "N vs. A"; and 83 overlapping DEGs were specifically found in "N vs. NA" and "A vs. NA". Therefore, a total of 666 unique up-regulated DEGs were identified and further selected for GO annotation and KEGG enrichment analysis (Figure 2b). Similarly, a total of 802 unique down-regulated DEGs were identified and further selected for GO annotation and KEGG enrichment analysis (Figure 2c).

Trend Co-Expression Analysis of DEGs
The 2319 DEGs were categorized into 15 gene clusters (Figure 3). A large proportion of the genes (1556, 67.1%) showed non-additive gene expression in livers of hybrid tilapia. Of these 1556 genes, 482, 112, 455, and 507 exhibited overdominance, high-parent dominance, low-parent dominance, and under-dominance patterns, respectively. Additionally, 763 DEGs exhibited additive-type gene expression, which accounted for 32.9% of all DEGs.

GO and KEGG Enrichment Analysis of DEGs
Based on GO annotation, 389 of 666 up-regulated DEGs were enriched into 42 GO terms (Supplementary Figure S2). Under biological process, the DEGs were enriched in those associated with metabolic processes and single biological processes. In cellular components, the DEGs were enriched in cellular components and cells. In molecular function, the DEGs were enriched in catalytic activity and binding.
By KEGG analysis, the up-regulated DEGs were significantly enriched in metabolic pathway, including fatty acid metabolism, carbon metabolism, peroxisome, biosynthesis of unsaturated fatty acids, fatty acid biosynthesis, glycine, serine and threonine metabolism, pyruvate metabolism, biosynthesis of amino acids, and arginine and proline metabolism (Figure 4a).  A total of 458 of the 802 down-regulated DEGs were enriched into 46 GO terms (Supplementary Figure S3). Under biological processes, the DEGs were enriched in cellular processes and metabolic processes. In cellular components, the DEGs were enriched in cellular fractions and cells. In molecular function, the DEGs were enriched in binding and catalytic activity. KEGG analysis showed that 802 down-regulated DEGs were significantly enriched in ribosome biogenesis in eukaryotes (Figure 4b).

Analysis of DEGs Potentially Related to Hybrid Tilapia Heterosis
Based on the GO and KEGG analysis, we selected 32 DEGs that may be most related to hybrid tilapia heterosis. The FPKM values of the DEGs for NA, N, and A are detailed in Table 2. We found that five genes related to fatty acid biosynthesis, including fatty acid synthase (FASN), long-chain fatty acid-CoA ligase (ACSL1, ACSL3, ACSL6), and acetyl coenzyme A carboxylase (ACACA) of NA had higher FPKMs value than both N and A, and so were long-chain fatty acid elongase (ELOVL5 and ELOVL6) and fatty acid desaturase (FADS2), involved in unsaturated fatty acid biosynthesis. The results provide further evidence that fatty acid metabolism enzyme genes were more active in hybrid tilapia.

Eukaryotic Biosynthesis and Defense-Related DEGs
The key enzymes of eukaryotic biosynthesis, SNU13, DKC1, Nug1/2, MPP10, UTP4, UTP5, UTP14, UTP10, and UTP15, were significantly down-regulated in hybrids compared with the parents. The antioxidant-related gene GSH-Px and HSP family genes (HSPA5 and HSP70) were also significantly down-regulated. Hierarchical cluster analysis of down-regulated DEGs is shown in Figure 5b.

Validation of RNA-Seq Data by qRT-PCR
As shown in Figure 6, the expression patterns of 16 DEGs detected by qRT-PCR were highly consistent with the data produced by RNA-Seq, which indicated that the sequencing results were reliable and accurate.

Discussion
Heterosis has been widely exploited in aquaculture breeding [3][4][5][6]. Hybrid tilapia, the offspring of maternal O. niloticus and paternal O. aureus, has been demonstrated to have enhanced performance traits, such as faster growth rate, and better stress and disease tolerance [14,15]. However, the molecular mechanisms underlying this phenomenon remain elusive. Zhou et al. [14] and Zhong et al. [15] illustrated the concept of GH and IGF-1 expression superiority and its mechanism for promoting vigorous growth in hybrid tilapia, respectively. In the current study, we constructed RNA-seq libraries of the livers of hybrid tilapia and its parents, and detected some DEGs and pathways that may be responsible for heterosis. Subsequently, trend co-expression analysis showed that non-additive gene expression patterns were prevalent in the liver of hybrid tilapia. In the end, based on the function of these DEGs and their location in the pathway, an overview of the superior growth generation was produced (Figure 7).

Analysis of Candidate DEGs Involved in Fatty Acid Metabolism
The levels of expression of numerous genes involved in the fatty acid metabolism pathway, such as FASN, ELOVL5, ELOVL6, ACSL6, ACSL1, ACSL3, and FADS2, were upregulated in the livers of hybrid tilapia compared with the parental species. Fatty acid synthase (FASN) is an important rate-limiting enzyme that catalyzes the synthesis of fatty acids used for energy generation, extended with C2 units derived from the elongation substrate malonyl-CoA [37]. The biosynthesis of long-chain polyunsaturated fatty acids requires the concerted activities of fatty acyl desaturase (FADS) and elongase (ELOVL) [38]. Similarly, Zheng et al. [18] found significant up-regulation of FAS, ELOVL1, ELOVL5, and ELOVL6 expression in a new backcross of T. gigas ♀ × (T. gigas ♀ × Bleak ♂) ♂, the growth advantage of this backcross may be related to the increased fatty acid synthesis. These up-regulated genes may lead to accelerated fatty acid synthesis and fatty acid metabolism.

Analysis of Candidate DEGs Involved in Carbon Metabolism
Carbon metabolism, which includes the EMP, PPP, and the TCA cycles, is very important for providing materials and energy for fish growth. Glycolysis plays a central role in the anabolism and catabolism of organisms [39], and provides intermediates for other metabolic pathways, such as gluconeogenesis, fatty acid metabolism, biosynthesis of amino acids, and the TCA cycle [40,41]. Almost all glycolytic genes were up-regulated in hybrid grouper compared with their parents [11]. In this study, both PFK1 and ENO1 in the glycolysis pathway were up-regulated in hybrid tilapia. The TCA cycle is a crucial pathway for generating energy through oxidizing carbohydrates and fatty acids [42]. In the current study, significantly enriched DEGs encoding ME3 and ACS were found in the TCA cycle. Additionally, the amino acid metabolic pathway genes GATM, SHMT, PSPH, and PSAT1, G6PD involved in PPP, were up-regulated in the F1 hybrids. These findings indicated a dramatic change in efficiency of amino acid metabolism between F1 hybrids and their parents. Consistent with the results of previous studies [8,31,43], we suggest that the enriched metabolism-related pathways identified in this study may directly contribute to growth heterosis in hybrid tilapia.

Analysis of Candidate DEGs Involved in Ribosome Biosynthesis in Eukaryotes
Ribosome biogenesis and protein synthesis are fundamental rate-limiting steps for cell growth and proliferation [44]. In eukaryotes, ribosome assembly is a rate-limiting step in ribosomal biogenesis. Ribosomal biogenesis is an elaborate process that determines the rate of protein synthesis, and by controlling a main anabolic cellular process, it regulates accumulation of cellular mass or cell growth [45]. Here, many key enzymes for eukaryotic ribosome biosynthesis in hybrid tilapia, SNU13, DKC1, Nug1/2, MPP10, UTP4, UTP5, UTP14, UTP10, and UTP15, were significantly decreased compared with their parents, which is consistent with the findings of Xiao et al. [7] on the domestic silkworm. Xie et al. [42] showed that the fast-growth phenotype of the ark shell Scapharca kagoshimensis may be attributed to lower energy requirements for metabolism maintenance, and more efficient protein biosynthesis and degradation. These findings revealed that efficiency of protein metabolism plays a role in growth heterosis in hybrid tilapia.

Analysis of Candidate DEGs Involved in Basal Defense Response
GSH-Px is an important peroxidase that belongs to the antioxidant defense system, which plays a fundamental role in overall defense mechanisms and strategies of organisms [46]. It was speculated that individuals with lower GSH-Px activity are susceptible to antioxidant protection [47]. HSPs are highly conserved molecular chaperones associated with helping other proteins fold properly during synthesis or repair of misassembled proteins [48]. HSPs, some defense genes and immune genes were significantly down-regulated in the largest individual sea cucumbers compared with the smallest ones [20]. In the current study, GSH-Px and HSPs (HSPA5 and HSP70) were significantly down-regulated in hybrid tilapia compared with their parents, indicating that its basal defense response was repressed. Pseudomonas syringae infection of Arabidopsis F1 hybrids demonstrated that the reductions in basal defense gene activity in these hybrids does not necessarily compromise their ability to mount a defense response comparable to the parents [49]. Consistent with the results of previous studies [50,51], the repressed expression of basal defense response may promote growth in hybrid tilapia.

Conclusions
Hybrid tilapia exhibited higher growth rates than their parents, O. niloticus and O. aureus. Here, we provide a global view of the transcriptomic divergence of hybrid tilapia and their parental lines using RNA-Seq. A total of 2319 DEGs were identified. Trend coexpression analysis revealed that non-additive expression was prominent in hybrid tilapia compared with their parents. Interestingly, the metabolic pathway was more highly regulated in fatty acid metabolism and carbon metabolism, whereas ribosome biosynthesis in eukaryotes and basal defense response were significantly down-regulated. These findings provide new insights into our understanding of growth heterosis in hybrid tilapia.

Supplementary Materials:
The following are available online at www.mdpi.com/article/10.3390/fishes7010043/s1. Figure S1: DEGs analysis and volcano plot for "NA vs. N", "NA vs. A", and "N vs. A". The x-axis is the value of Log2 (Fold Change), and the y-axis is the value of -Log10 (FDR). The red dots reveal the up-regulated DEGs, the green dots reveal the down-regulated DEGs, and black dots reveal non-differentially expressed genes. NA, N, and A denote hybrid tilapia, O. niloticus, and O. aureus, respectively. Figure S2: GO annotation of significantly up-regulated DEGs between O. niloticus or O. aureus and hybrid tilapia. Figure S3: GO annotation of significantly down-regulated DEGs between O. niloticus or O. aureus and hybrid tilapia. Table S1. Summary statistics of clean transcriptome sequencing data. Table S2. Sequencing data and selected reference genome sequence alignment results.