Integrated Analysis of Transcriptomic, miRNA and Proteomic Changes of a Novel Hybrid Yellow Catfish Uncovers Key Roles for miRNAs in Heterosis*

mRNA-seq, miRNA-seq, proteomes of P. fulvidraco, P. vachelli, the hybrid Huangyou-1 livers were characterized. The nonadditive, HEB and ELD pattern were readily identified in transcriptional, post-transcriptional, or protein levels. Several predicted miRNA-mRNA-protein pairs were found and validated by qRT-PCR and PRM. Immune defense, metabolism, digestion and absorption, cell proliferation and development involved in the generation of the heterosis phenotype. The high parental genes/proteins coupled with low parental miRNAs of the offspring are inherited from the mother or father. Graphical Abstract Highlights mRNA-seq, miRNA-seq, proteomes of P. fulvidraco, P. vachelli, hybrid Huangyou-1. Predicted miRNA-mRNA-protein pairs were found and validated by qRT-PCR and PRM. Immune, metabolism, digestion, absorption, proliferation, development generate heterosis. High parental gene/protein with low parental miRNAs inherit from the mother or father. Heterosis is a complex biological phenomenon in which hybridization produces offspring that exhibit superior phenotypic characteristics compared with the parents. Heterosis is widely utilized in agriculture, for example in fish farming; however, its underlying molecular basis remains elusive. To gain a comprehensive and unbiased molecular understanding of fish heterosis, we analyzed the mRNA, miRNA, and proteomes of the livers of three catfish species, Pelteobagrus fulvidraco, P. vachelli, and their hybrid, the hybrid yellow catfish “Huangyou-1” (P. fulvidraco ♀ × P. vachelli ♂). Using next-generation sequencing and mass spectrometry, we show that the nonadditive, homoeolog expression bias and expression level dominance pattern were readily identified at the transcriptional, post-transcriptional, or protein levels, providing the evidence for the widespread presence of dominant models during hybridization. A number of predicted miRNA-mRNA-protein pairs were found and validated by qRT-PCR and PRM assays. Furthermore, several diverse key pathways were identified, including immune defense, metabolism, digestion and absorption, and cell proliferation and development, suggesting the vital mechanisms involved in the generation of the heterosis phenotype in progenies. We propose that the high parental expression of genes/proteins (growth, nutrition, feeding, and disease resistance) coupled with low parental miRNAs of the offspring, are inherited from the mother or father, thus indicating that the offspring were enriched with the advantages of the father or mother. We provide new and important information about the molecular mechanisms of heterosis, which represents a significant step toward a more complete elucidation of this phenomenon.

Heterosis is a complex biological phenomenon in which hybridization exhibits superior phenotypic characteristics such as enhanced growth rate, development, and stress tolerance relative to the parents (1). Dominance, overdominance, and epistasis have each been proposed as explanations of heterosis from a genetic perspective. However, these hypotheses are largely conceptual and are unrelated to molecular principles; thus, they fall short in explaining the molecular basis of heterosis (2). Although heterosis has been exploited extensively in fish production for many decades, the molecular mechanisms underlying the phenomenon remain largely unknown despite more than a century of study (3)(4)(5). The utilization rate of fish heterosis has far exceeded understanding of it on a theoretical level. This lack of understanding also limits improvements to its use in aquaculture.
Hybridization is a fundamental process in evolution that results in the emergence of novel genotypes from the merging of two different genomes. It is theorized that because no new genes are produced in the hybrids, heterosis is likely caused by differences resulting from qualitative or quantitative modifications of gene expression (6 -8). Several studies of plants using quantitative real-time PCR (qRT-PCR) 1 and next-generation sequencing technologies (NGST) have enabled detailed investigations of allelic heterozygosity and/or the epigenetic changes resulting in heterosis in hybridizations relative to noncoding RNA (9), methylation (10), and transcriptome changes (11). Recent studies of hybrid cyprinidae (Megalobrama amblycephala & ϫ Culter alburnus ( (12) and Carassius auratus & ϫ Cyprinus carpio ( (13)) have focused on the nonadditively expressed genes (NEG), homoeologue expression bias (HEB), or expression level dominance (ELD) patterns to summarize gene regulation and their underlying mechanisms analyzed by NGST, such as the superiority of cell development, stress, and adaptability. Also, several elevated gene expressions detected by qRT-PCR in hybrid fish have been reported, which may affect growth rate or disease resistance (14,15). MicroRNAs (miRNAs) play a critical role in the regulation of gene expression at the post-transcriptional level, especially in epigenetic gene regulation. Data from several plant studies have shown that epigenetic variation by cis-or trans-regulation of miRNAs may contribute to the molecular mechanisms of complex traits such as growth (9), vigor (16), stress response (17), and heterosis (18). The miRNA expression in hybrid cyprinidae leads to nonadditive expression of target genes that may affect cellular functional stability in ploidy change, stress adaption, and development (12,19). The genome-wide gene expression bias of miRNA in interspecific hybrid fishes is essential for heterosis. However, the effects of miRNA-mRNA pairs for hybrid vigor have yet to be illustrated. Integrated analysis using NGST has been performed on several teleosts for a better understanding of the mRNA-miRNA network of intermuscular bone development (20), the sex determination (21), hypoxia adaptation (22), and albinism (23). A potentially more reliable method for predicting the miRNA-mRNA relationship within a particular biological context is to integrate real mRNA and miRNA transcriptomic data into in silico target predictions.
The transcriptional or post-transcriptional level of research on the interpretation of the genetic mechanism of heterosis has improved significantly (24). Nevertheless, transcription change does not depict the various changes of protein levels with complete accuracy because of the limitations of transcription regulation (25). Therefore, from the level of proteomics study of the differences between hybrids and their parent proteins, the change rule for understanding differences in proteins in species plays a role in heterosis, where an in-depth discussion is indispensable. Recently, proteome studies of different plant organs using two-dimensional fluorescence difference gel electrophoresis have shed light on the potential molecular mechanisms of heterosis (26 -29). The most commonly observed patterns indicate the global expression trend of both nonadditivity and additivity in regulating plant heterosis (8,30). Thus far, no research has examined the proteomic profiles of fish heterosis.
In the past decades, various catfish species have been widely cultivated in many countries, and the industry has developed rapidly (31). To resolve issues of variable quality and the degradation of economic characteristics, hybridization technology has been introduced into the artificial breeding of catfish (32)(33)(34). It is worth noting that the artificial hybrid yellow catfish Huangyou-1 (hybrid F1, H), produced from the hybridization of Pelteobagrus fulvidraco & (female, F) ϫ P. vachelli ( (male, M), has advantages such as a seductive body color, faster growth and development, a high feeding rate, and greater immunity when compared with its parents (35). These special characteristics of H suggest that it is not only a significant aquaculture species but also a potential model organism for study of the molecular mechanisms of heterosis. At present, many research studies have examined the biological and physiological features of H (35)(36)(37); however, very little is known about the genetic mechanisms of heterosis.
The present NGST and Tandem Mass Tags (TMT) study should provide unprecedented resources to address such questions as to how hybridization affects gene/protein expression and changes the molecular pathways, which could lead to broader adaptability in nascent hybrid fish, and whether miRNAs participate in this process. For this purpose, we performed mRNA-seq, miRNA-seq, and proteomics in liver from H, F, and M ( Fig. 1), which is the first report on the integrated analysis of fishes and offers deeper insight into heterosis miRNA-mRNA-protein pairs, also validated by qRT-PCR and parallel reaction monitoring (PRM) assays (38). Our data suggest a potential miRNA moderated regulatory mechanism for gene expression as well as protein expression as the basis for heterosis phenotypes in hybridization.

EXPERIMENTAL PROCEDURES
Hybridization and Cultivation-Parental fish of P. fulvidraco and P. vachelli were obtained from Nanjing Fisheries Research Institute. Fish breeding and cultivation were performed at Nanjing Normal University in Jiangsu Province, China. Artificial fertilization of M & ϫ M (, F & ϫ F (, and F & ϫ M ( were conducted through manual abdominal compression for harvesting eggs from females and stripping sperm from males after mixed spawning-stimulating hormone injections (female parent: 2300 IU kg-1 human chorionic gonadotropin, 22 g/kg luteinizing hormone releasing hormone A2, 1.6 mg/kg domperidone, and 0.2 mg/kg carp pituitary; male parent: half dosage; Sansheng pharmaceutical Co. Ltd., Henan, China).
Obtained fertilized eggs (M, F, and H) were incubated at a water temperature of 28°C in flow-through water until hatching (about 2 days postfertilization). Three lines were fed with live Artemia nauplii to satiation thrice daily until the mouth gape sizes for all fish matched the size of the artificial compound feed, which composed of Ͼ 42.0% protein, 6.0% fiber, 14% ash, 0.5%ϳ3.5% calcium, 1.2% phosphorus and 0.3%ϳ2.5% sodium chloride (Jiaji Feed Co. Ltd., Zhenjiang, China). The feed was slightly oversupplied based on previous consumption. For M, 600 juvenile fish were averagely reared in 3 aquaria with the same stocking density (200 individuals/1350 L of volume with bio-filtered water recirculation system; flow rate of 45 L/min; same settings were applied to all aquaria). The treatment procedures for F and for H were the same with M. If there were dead fish in these aquaria, we replaced them with specific individuals of the same weight.
During the experiment period, all fish were supplied with tap water aerated for more than 24 h; the photoperiod was maintained at 11 h light/13 h darkness; the water temperature was monitored daily at 8:00; and other water quality characteristics were monitored weekly. Water quality variables of all aquaria were as follows: total ammonia (NH 4 ϩ ϩNH 3 ) of 0.2-1.0 mg/L; nitrite of 0.1-0.5 mg/L; pH of 6.6 -7.4; dissolved oxygen Ͼ 4.0 mg/L. Experimental Design and Statistical Rationale-After cultivation at 28 Ϯ 1°C for 105 days and sequential feed restriction for 2 days, the total body weight of 200 juveniles in each aquaria (total 9 aquaria) were monitored. Mid-parent heterosis (MPH, %) ϭ 100 ϫ (F1-Pn)/Pn, where F1 is the mean performance of the hybrids and Pn is the mean value of both parents. For exploring heterosis, 6 individuals were randomly chosen from each aquaria (total 9 aquaria). Samplings of M, F, and H had three biological replicates (M1, M2, M3, F1, F2, F3, H1, H2, and H3) from the three aquaria, with each sample comprised of a pool of above 6 different individual liver tissues. Fish were killed by dissection of livers after mild anaesthetization in a eugenol bath (1:10,000). These harvested samples were used to conduct mRNA-seq, miRNA-seq, proteomics, qRT-PCR, and PRM, respectively (Fig. 1). All the above experiments including the mass spectrometry experiments had three biological replicates (M1, M2, M3, F1, F2, F3, H1, H2, and H3).
Statistical analysis of qRT-PCR and PRM (Tier 2 level) (38) was performed using SPSS 22.0. The relevant values of three groups were analyzed by one-way analysis of variance test, followed by an unpaired, two tailed t test. Statistical significance was considered at p value Ͻ 0.05 (39). All data were expressed as mean Ϯ standard errors in terms of relative expression.
Differentially Expressed Genes (DEGs) Analysis-Salmon (40) was used to perform expression level for unigenes by calculating Transcripts Per Million (TPM) (41). The DEGs were selected with a fold change Ͼ 2 and with statistical significance (FDR Ͻ 0.01) by R package edgeR (42,43).
Analysis of Differentially Expressed miRNAs (DEMs)-miRNA differential expression based on normalized deep-sequencing counts was analyzed by selectively using a student's t test based on the experimental design (22). The significance threshold in this test was set to be p value (P) Ͻ 0.05 (21,44).
The Quantitative Data Analysis of Proteomics-For multiple comparisons of statistical analysis between M and F, nine protein ratios could be produced after comparisons among three M and three F (i.e. M1/F1, M2/F1, M3/F1, M1/F2, M2/F2, M3/F2, M1/F3, M2/F3, and M3/F3). After this, the p values were generated via permutation test. Fold change was calculated as the mean comparison of pairs among nine protein ratios (45). Accordingly, we set the ratio threshold at Ͼ1.2 for protein up-regulation and at Ͻ0.83 (1/1.2) for protein downregulation (46). Proteins with a fold change larger than 1.2 and a p value Ͻ 0.05 were significantly differentially expressed (47,48) Transcriptome Sequencing-Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's protocol. The total RNA quantity and purity were analyzed with Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA) with RIN number Ͼ7.0.
For nine cDNA library constructions (M1, M2, M3, F1, F2, F3, H1, H2, and H3), ϳ10 g of total RNA representing a specific adipose type was subjected to isolate Poly (A) mRNA with poly-T oligo attached magnetic beads (Invitrogen). Following purification, the mRNA was fragmented into small pieces using divalent cations under elevated temperature. Then, the cleaved RNA fragments were reversetranscribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, CA), the average insert size for the paired-end libraries was 300 bp (Ϯ50 bp). Paired-end sequencing was then performed on an Illumina Hiseq 2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol (49).
De Novo Assembly, Unigene Annotation and Functional Classification-First, Cutadapt (50) and perl scripts in house were used to remove the reads that contained adaptor contamination, low quality bases, and undetermined bases. Then sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) and included the Q20, Q30, and GC-content of the clean data. All downstream analyses were based on clean data of high quality. De novo assembly of the transcriptome was performed with Trinity 2.4.0 (51). Trinity groups were transcripted into clusters based on shared sequence content. Such a transcript cluster is very loosely referred to as a "gene." The longest transcript in the cluster was chosen as the "gene" sequence (aka Unigene).
Small RNAs Sequencing-Approximately 5 g of total RNA were used to prepare nine small RNA libraries (M1, M2, M3, F1, F2, F3, H1, H2, and H3) according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA), and then the libraries were sequenced by Illumina Hiseq 2500 at the LC-BIO following the vendor's recommended protocol. Raw reads were subjected to an in-house program, ACGT101-miR (LC Sciences, Houston, TX), to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeats.
Subsequently, unique sequences with lengths of 18ϳ26 nucleotides were mapped to specific species precursors in miRBase 21.0 by BLAST search to identify known miRNAs and novel 3p-and 5pderived miRNAs. Length variation at both 3Ј and 5Ј ends and one mismatch inside of the sequence were allowed in the alignment. The unique sequences mapped to specific species with mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapped to the other arm of known specific species with precursor hairpins opposite the annotated mature miRNA-containing arm were novel 5p-or 3p-derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 21.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the specific species genomes to determine their genomic locations. The above two we defined are known as miRNAs. The unmapped sequences were BLASTed against the specific genomes, and the hairpin RNA structures containing sequences were predicated from the flank 80 nt sequences using RNAfold software (http://rna.tbi. univie.ac.at/cgi-bin/RNAfold.cgi). The criteria for secondary structure prediction were the following: (1) the number of nucleotides in one bulge in stem (Յ12), (2) the number of base pairs in the stem region of the predicted hairpin (Ն16), (3) the cutoff of free energy (kCal/mol ՅϪ15), (4) the length of the hairpin (up and down stems ϩ terminal loop Ն50), (5) the length of the hairpin loop (Յ20), (6) the number of nucleotides in one bulge in the mature region (Յ8), (7) the number of biased errors in one bulge in the mature region (Յ4), (8) the number of biased bulges in the mature region (Յ2), (9) the number of errors in the mature region (Յ7), (10) the number of base pairs in the mature region of the predicted hairpin (Ն12), and (11) the percentage of mature region in the stem (Ն80).
Prediction of Target Genes of miRNAs-To predict the genes targeted by DEMs, two computational target prediction algorithms (Tar-getScan 50 and miRanda 3.3a) were used to identify miRNA binding sites. Finally, the data predicted by both algorithms were combined and the overlaps were calculated.

mRNA, miRNA and Protein Changes in Hybrid Catfish Heterosis
TMT Proteomic Analysis-Liver samples (H1, H2, H3, F1, F2, F3, M1, M2, and M3) were first frozen to a dry powder with a vacuum freeze drier. The freeze-dried powder was dissolved with 200 l TEAB (Triethylammonium bicarbonate) dissolution buffer (8 M Urea/100 mM TEAB), and then after centrifugation at 16,000 ϫ g (the centrifugal radius was 0.1 m) for 20 min, the supernatant subsided by adding 4-fold volume cold acetone containing 10 mM dithiothreitol (DTT) for about 2 h. The resuspended powder was incubated at Ϫ20°C for ϳ2 h. After centrifugation at 16,000 ϫ g for 20 min at 15°C, the precipitate was collected and mixed with 800 l cold acetone containing a solution of 10 mM DTT for 1 h at 56°C to break the proteins' disulfide bonds. Following a second round of centrifugation, this time at 16,000 ϫ g for 20 min at 15°C, the precipitate was collected, dried, and then stored at Ϫ80°C for later use.
Total protein concentration was measured using the Bradford method (53). For each sample, 100 g of protein was dissolved to 500 l in a dissolution buffer, and then diluted with 500 l 50 mM NH 4 HCO 3 . The samples were reduced with 10 mM DTT at 56°C for 30 min, alkylated with 50 mM iodoacetamide for 30 min in the dark, and diluted 4 times with 10 mM TEAB.
The peptide mixture was redissovled in the buffer A (buffer A: 10 mM ammonium formate in water, pH 10.0, adjusted with ammonium hydroxide), and then fractionated by high pH separation using an Aquity UPLC system (Waters Corporation, Milford, MA) connected to a reverse phase column (BEH C18 column, 2.1 mm ϫ 150 mm, 1.7 m, 300 Å Waters Corporation, Milford, MA). High pH separation was performed using a linear gradient. Starting from 0% B to 45% B in 35 min (B: 10 mM ammonium formate in 90% ACN, pH 10.0, adjusted with ammonium hydroxide). The column flow rate was maintained at 250 l/min and column temperature was maintained at 45°C. Twelve fractions were collected, and each fraction was dried in a vacuum concentrator for the next step.
Low pH nano-HPLC-MS/MS Analysis-The fractions were resuspended with 32 l solvent C respectively (C: water with 0.1% formic acid; D: ACN with 0.1% formic acid), separated by nanoLC and analyzed by online electrospray tandem mass spectrometry. The experiments were performed on an EASY-nLC 1000 system (Thermo Fisher Scientific, Waltham, MA) connected to an Orbitrap Fusion mass spectrometer (Thermo Fisher Scientific, San Jose, CA) equipped with an online nano-electrospray ion source. A 4 l peptide sample was loaded onto the trap column (Thermo Scientific Acclaim PepMap C18, 100 m ϫ 2 cm), with a flow of 10 l/min for 3 min and subsequently separated on the analytical column (Acclaim PepMap C18, 75 m ϫ 25 cm) with a linear gradient, from 5% D to 30% D in 95 min. The column was re-equilibrated at initial conditions for 15 min. The column flow rate was maintained at 300 nL/min and column temperature was maintained at 45°C. The electrospray voltage of 2.0 kV versus the inlet of the mass spectrometer was used.
The Orbitrap Fusion mass spectrometer was operated in the datadependent mode to switch automatically between MS and MS/MS acquisition (56). Survey full-scan MS spectra (m/z 400 -1,600) were acquired in the Orbitrap with a mass resolution of 60,000 at m/z 200. The automatic gain control (AGC) target was set to 500,000, and the maximum injection time was 50 ms. MS/MS acquisition was performed in the Orbitrap with a 3 s cycle time, and the resolution was 15,000 at m/z 200. The intensity threshold was 50,000, and the maximum injection time was 150 ms. The AGC target was set to 150,000 and the isolation window was 2 m/z (57,58). Ions with charge states 2ϩ, 3ϩ, and 4ϩ were sequentially fragmented by higher energy collisional dissociation (HCD) with a normalized collision energy (NCE) of 37%. In all cases, one microscan was recorded using dynamic exclusion of 30 s. MS/MS fixed first mass was set at 110.
The Database Searching-The original mass spectrum data was converted into MGF format by the integration tool MSConvert of ProteoWizard (Version 2.1). All MS/MS samples were analyzed using Mascot (Matrix Science, London, UK; version 2.3) and Percolator software. The Mascot was set up to search the transcriptome data in this paper, and the RNA-seq results were employed for database searching (Accession codes: GSE118627; entries: 17,733). Mascot was searched with a fragment ion mass tolerance of 0.050 Da and a parent ion tolerance of 10 PPM, with allowance for one missed cleavage upon trypsin digest. Several parameters in Mascot were set for peptide searching, including oxidation of methionine as potential variable modifications and carbamidomethyl (C), TMT 10-plex (Nterm) and TMT 10-plex (K) as fixed modifications. Use the percolator algorithm to control peptide level false discovery rates (FDR) lower than 1% (45). Only unique peptides were contained for TMT labeling quantification, and the method of normalization on protein median was used to correct experimental bias, the minimum number of proteins that must be observed to allow was set to 1000. Several parameters in IQuant software were set for protein quantification, including all unique peptide as quant peptide, at least one unique spectrum as quant number, variance stabilization normalization as normalization standard, weighted average as protein ratio, permutation test as statistical analysis technique (59).
Analyses of Nonadditivity/additivity, Homoeolog Expression Bias (HEB), and Expression Level Dominance (ELD)-Values from the three biological replicates of each parent were averaged to calculate the MPV, and then the expression levels of F1 significantly (p Ͻ 0.05, Fisher's exact tests) different from the MPV were designated as nonadditive.
To identify putative orthologues between F and M, the two sets of sequences were aligned using the reciprocal BLAST (BLASTN) hit method, with an e-value cutoff of 1e Ϫ20 (60). Two sequences were defined as orthologues if each were the best hit of the other and if the sequences were aligned over 300 bp. After identifying SNPs between the F and M orthologs, we mapped our reads from F and M to compare the mapping results. Reads with SNPs that differed between the F-and M-genome in the progenitors were parsed into F and M homoeolog-specific bins using custom perl scripts (15). For HEB analyses, we compared F and M homoeolog expression levels in the three hybrids. The fold change Ͼ 2 and with statistical significance (FDR Ͻ 0.01) of homoeolog expression value had been considered to be the threshold value of F-HEB. In contrast, the FC Ͼ 2 and FDR Ͻ 0.01 had been considered the threshold values of M-HEB (13).
In the hybrid offspring, genes/miRNAs/proteins that were identified as differentially expressed in the hybrid relative to the diploid parents were binned into 12 possible differential expression categories as described in Shen et al. (2017) (18), with the term "additivity" changed to "Middle parents (additivity or partial-dominance)." We further designated genes/miRNAs/proteins with expression levels in F1 that were statistically like those in the male and female parent as ELD-M and ELD-F genes/miRNAs/proteins, respectively.
KEGG Pathway Enrichment Analysis-To identify candidate biomarkers, we employed the hypergeometric test to perform KEGG pathway enrichment (44). KEGG pathway enrichment analysis was used to identify significantly enriched KEGG pathways using the corrected p value Ͻ 0.05 as a threshold for finding significantly enriched KEGG terms in the input list of DEGs and DEPs, comparing them to the whole genome background. The calculation formula of p value was as follows: N represented the number of KEGG annotated proteins in P. vachelli, n represented the number of DEGs or DEPs in N, M represented the number of particular KEGG annotated genes or proteins in a genome, and m represented the number of particular KEGG annotated genes or proteins expressed differentially in M. In this article, all graph data for KEGG pathway were significantly enriched data.
Integrated Analysis of Transcriptome, miRNAome, and Proteome-To define all the possible miRNA-mRNA interactions, including positive and negative relationships between miRNA and mRNA expression, we used ACGT101-CORR 1.1 to construct the miRNA-mRNA regulatory network (22). Briefly, we normalized all the sample-matched miRNA and mRNA sequencing data. Afterward, integration of miRNA-seq with mRNA-seq was achieved by integrating expression profiles of DEGs and DEMs with the addition of differentially expressed miRNA-targeting information (21). Because protein identification was employed against the above RNA-seq results for database searching, the proteins added miRNA-mRNA pairs directly to the miRNA-mRNA-protein.
Predicted miRNA-mRNA-protein Validation by qRT-PCR-The expression profiles of 6 DEM (ssa-miR-30e-5p, bta-miR-221, dre-miR-130c-5p, pma-miR-181a-3p, chi-miR-140 -3p, and cgr-miR-27a-3p_Rϩ1_1ss19CT), 6 DEG and 6 DEP among the miRNA-mRNA-protein interaction networks were further validated using qRT-PCR and PRM (t test, p Ͻ 0.05), bile salt export pump (Abcb11), Putative all-trans-retinol 13,14-reductase (Retsat), Calcium/calmodulin-dependent protein kinase type II delta 1 chain (Camk2d1), Truncated complement factor I (Cfi), Neural cell adhesion molecule L1 (CHL1), and acyl-coenzyme A thioesterase 2, mitochondrial-like (ACOT6). The relative expression of 6 mature miRNAs were selected and analyzed by quantifying the miRNA stem-loop. Total RNAs were isolated using Trizol reagent (Biotake, Beijing, China) following the cDNA generation using 1 g of total RNA by reverse transcription kit (Toyobo, Osaka, Japan). Quantitative real-time PCR (qRT-PCR) was performed on an ABI Step One Plus system (Applied Biosystems, Foster, CA) by using qRT-PCR reagents provided by Toyobo. The stem-loop primers are shown in the supplementary data (supplemental Table S1). PCR amplification was conducted under an initial denaturation at 94°C for 30 s, and then 40 cycles of amplification including the denaturation at 94°C for 20 s, annealing at 61°C for 30 s, and extension at 72°C for 30 s; after 40 cycles was the final extension at 72°C for 1 min. Each sample was tested in triplicate. Finally, the melting curve was performed to verify the specificity of PCR amplification. U6 was used as an internal control (21,22).
To verify RNA-seq results, qRT-PCR method with ␤-actin as an internal control was used to explore mRNA expression levels. qRT-PCR was performed with an SYBR Green Master kit according to the manufacturer's protocol (Roche, Basel, Switzerland) (45). The primers for qRT-PCR are listed in supplemental Table S1. The experiments were carried out in triplicate with a total volume of 20 l using ABI stepone TM plus (Applied Biosystems, Waltham, MA), containing 10 l of SYBR green master, 4 l of cDNA (500 ng), and 3 l of forward and reverse primers (2 mol/L). The qRT-PCR was programmed at 95°C for 10 m, followed by 40 cycles of 95°C for 15 s, and 55°C for 1 min. The expression level was calculated by 2 Ϫ‚‚CT method and subjected to statistical analysis. During mRNA expression analysis, foldchange was determined by comparing with the reference gene expression.

Development and Analytical Validation Targeted MS Assays-Mul-
tiple reaction monitoring (MRM) is broadly applicable for verification phase Tier 2 multiplexed multiple reaction monitoring assay development within the Food and Drug Administration-National Cancer Institute biomarker development pipeline (61,62). Parallel reaction monitoring (PRM) provides a mass analyzer with higher resolution, better quality, and greater precision than MRM (39,63). To further verify the protein expression levels obtained by TMT free analysis, the expression levels of selected proteins were further quantified by LC-PRM-MS analysis (Tier 2 level) in Shanghai Applied Protein Technology Co., Ltd. These proteins were manually selected as representatives for their potential roles in heterosis. Briefly, peptides were prepared according to the TMT free protocol, and an absolute quantification stable isotope peptide was spiked in each sample as internal standard reference. As internal standards, we synthesized DSP-SAPVNVTVR stable isotope-labeled peptides, which included a valine enriched in stable isotopes (six 13C atoms). Tryptic peptides were loaded on C18 stage tips for desalting prior to reversed-phase chromatography on an Easy nLC-1200 system (Thermo Scientific). Onehour liquid chromatography gradients with acetonitrile ranging from 5 to 35% in 45 min were used (38).
PRM analysis was performed on a Q Exactive Plus mass spectrometer (Thermo Scientific) (64). Methods optimized for collision energy, charge state, and retention times for the most significantly regulated peptides were generated experimentally using unique peptides of high intensity and confidence for each target protein. The targeted peptides were isolated with a 2Th window. Ion activation/dissociation was performed at normalized collision energy of 27 in a higher energy dissociation (HCD) collision cell. The raw data were analyzed using Skyline (MacCoss Lab, University of Washington) where signal intensities for individual peptide sequences for each of the significantly altered proteins were quantified relative to each sample and normalized to standard reference. Three sub-ions with high abundance and as continuous as possible were selected for quantitative analysis, and the peak area of each target peptide segment was derived after Skyline analysis (65).
Ethics approval and consent to participate-All experiments were performed according to the Guidelines for the Care and Use of Laboratory Animals in China. This study was also approved by the Ethics Committee of Experimental Animals at Nanjing Normal University for Grant No. SYXK (Jiangsu) 2015-0028.

Mid-parent Heterosis (MPH, %) of Hybrid Yellow Catfish
Huangyou-1 (H)-In the process of cultivation at 28°C, H showed obvious heterosis at about 90 days. After cultivation at 28 Ϯ 1°C for 105 days with the same environment and sequential feeding restrictions for 2 days, the body weight of M, F, and H were 10.83 Ϯ 0.82, 9.44 Ϯ 0.52, and 13.11 Ϯ 0.98 g and the body length of M, F, and H were 8.60 Ϯ 0.42, 7.78 Ϯ 0.41, and 9.51 Ϯ 0.53 cm (with the standard error the mean of the three aquaria). Mid-parent heterosis (MPH, %) for body weight and body length were 29.35% and 16.12%, respectively. These heterosis results were like findings of previous research (36,37). We used liver tissues from three lines (M, F, and H) to investigate the po-tential molecular mechanisms underlying H heterosis by NGST and TMT analysis (Fig. 1).
Summary of Transcriptome, miRNAome, and Proteome-All MS/MS spectra were processed using ProteoWizard, Mascot and Percolator software. TMT analysis of P. vachelli liver proteome showed 544,168 spectra in this database (75,136 identified spectra) and resulted in 29,562 identified peptides. In summary, a total of 57,710 unique genes, 818 unique miRNAs (supplemental Table S2), and 4528 unique proteins had been identified (Fig. 2); quality control analysis of the transcriptome, miRNAome, and proteome is shown in supplemental Table S3, including TMT labeling efficiency of peptide segment, The distribution diagram of missed cleavages, Molecular weight profile of protein, The distribution diagram of Peptide length, Specific peptide number distribution map (supplemental Table S3).
The hierarchical clusterings of the DEP, DEM, and DEG based on the fold change' Z score of M1/(M2, M3, F1, F2, F3, H1, H2, H3), nine samples' log10 (Norm) and nine samples' Z score of TPM were showed in supplemental Fig. S1. The results showed that the nine samples were sorted into three distinct groups based on three species (supplemental Fig.  S1). The peptide/protein identifications and quantitation for all proteins were shown in supplemental Table S4. All coefficient of variation (CV) distribution for proteome showed that mean CV of F versus H, M versus H, M versus F were 0.07, 0.062 and 0.065, respectively, and proteome project has a good repeatability (supplemental Fig. S2).
Crucial Additive Genetic Pattern for H Heterosis-Dynamic changes of the transcriptome, miRNAome, and proteome in F1 and its two parents were investigated by analyzing the liver of three lines (M, H, and F) cultivated at 28°C for 105 days (Fig. 1). Genes, miRNAs, and proteins in F1 that had expression levels significantly different from the mid-parent value (MPV; the average expression value of the two parents; p Ͻ 0.05, Fisher's exact test) were designated as NEG, nonadditively expressed miRNAs (NEM), and nonadditively expressed proteins (NEP) (66). We found that a small portion (Ͻ9%) of all expressed genes and proteins were nonadditive, with miRNAs having a slightly greater portion (22.5%). This suggests that the majority of genes, miRNAs, and proteins in H show additive expression. This crucial additive genetic pattern provided an insight into the character of the hybrid related to heterozygosity, in which two different alleles from different species cooperate in the control of regulatory function. Only a small fraction of the expressed genes or proteins showed nonadditivity of expression (Fig 2A, supplemental Table S5), and such nonadditive genes, miRNAs, and proteins could potentially be responsible for generating the heterosis phenotype.
By performing KEGG pathway analyses for NEG, a total of 40 pathways were significantly (p Ͻ 0.05) identified (Fig. 2B). Among these pathways, the functions of immune defense, antioxidant capacity, metabolism (vitamins, amino acids, lipids, and carbohydrates), and digestion and absorption were the most frequently enriched. For NEP, 12 pathways, including digestion and absorption and immune defense, were enriched (Fig. 2B). Although some pathways related to plant and human disease were enriched by H NEG/NEP, it implies that the functions (e.g. disease resistance) of these pathways are included in H hybrid vigor.

FIG. 1. Work flow of bioinformatics analysis conducted to investigate the heterosis mechanism in hybrid yellow catfish Huangyou-1 (P. fulvidraco & ؋ P. vachelli () in terms of three main aspects (i.e. changes in transcriptional, post-transcriptional, and protein levels).
Abbreviations: DEG, differentially expressed genes; DEM, differentially expressed miRNAs; DEP, differentially expressed proteins. Fig. 3 shows summaries of differentially expressed genes (DEG), differentially expressed miRNAs (DEM), and differentially expressed proteins (DEP) in all comparisons. The percentage of genes/miRNAs/proteins showing differential expressions between F1 and their two parents were asymmetric (p Ͻ 0.05, Fisher's exact test). DEG, DEM, and DEP between the M and F lines were obviously higher than those of M-H and F-H. For example, ϳ12,786 of genes were differentially expressed between M and F, whereas the DEG numbers of M-H and F-H were 6808 and 4593 (Fig. 3A).

DEG, DEM, and DEP in the Two Parents Enriched in the Hybrid F1-
Interestingly, we found that a high percentage DEG/DEM/ DGP of the M-H and F-H overlapped with M-F (Fig. 3). Genes, miRNAs, and proteins differing in the expression level between the two parents were overrepresented (Ͼ67%) among the genes differentially expressed in the hybrids and parents, suggesting that the DEG/DEM/DEP of M-F were enriched in those of M-H or F-H in transcriptional, post-transcriptional, and protein levels. Notably, this phenomenon was more obvious (Ͼ90%) at the transcriptional level (Fig. 3A). With all of the above results, the differences between offspring and parents rooted in the differences among parents suggested a genetic pattern that tends to parental HEB and ELD in H heterosis.
Vital Genetic Pattern of Homoeolog Expression Bias (HEB) in Hybrid F1-To address whether the observed category of HEB truly reflects the homoeolog expression bias in the hybrid yellow catfish Huangyou-1, we compared the parental diploids and diploid hybrids of 11,789 homologous genes with homoeolog-specific SNPs on a case-by-case basis (supplemental Table S6).
For expression in parents, Rows 1, 6, and 7 (77.8%) in Table  I show that an equal pattern of parental expression was maintained for more than half of all homoeologous genes in this analysis (p Ͻ 0.05 in comparisons, Fisher's exact test). For expression in progeny, many genes (Rows 6, 7, 8, and 9; p Ͻ 0.05 in comparisons, Fisher's exact test) were detected as having novel patterns that accompanied the genome merger. These cases suggested novel regulatory and/or evolutionary interactions in the hybrid offspring. Rows 1, 4, and 5 represent the second most common class of genes, the parental expression patterns, representing 13% of the 11,789 genes. Rows 2 and 3 (1.6%) represent a small number of the pre-existing expression bias in the parental homoeolog reverted to non-differential expression of the homoeologous copies in the diploid hybrids. Significant HEB in hybrid yellow catfish Huangyou-1 (Rows 11 and 12) showed the expression patterns of homologous genes in progeny that favor the paternal fish. We also collected genes with overall no or F/M bias in progeny in hybrid yellow catfish Huangyou-1 (Rows 13 and 14). Our result for H showed 95.4% (Table I,

mRNA, miRNA and Protein Changes in Hybrid Catfish Heterosis
Importance of Parental ELD For Heterosis-To dissect the differential expression pattern in the hybrid F1 (H) relative to the parents (F or M), we classified genes into 12 male-hybridfemale (M-H-F) expression patterns as described in Shen et al. (18,67,68), with the term "additivity" changed to "Middle parents (additivity or partial-dominance)" (Fig. 4A). Interestingly, many parental ELD patterns were identified, with more than 80% (7093/8833) in transcriptome DEG, 70% (85/121) in miRNAome DEM, and 65% (785/1,200) in proteome DEP, respectively (Fig. 4A, supplemental Table S7-S9). These contrasted with the limited number of genes found in these samples that exhibited middle parents and transgressive regulation. The number of ELD-F pattern (types III and IV) was significantly higher than that of the ELD-M pattern (types V and VI) in transcriptome DEG, miRNAome DEM, and proteome DEP (p Ͻ 0.05, Fisher's exact test). These findings indicate that parental ELD patterns play an important role in fish allodiploid of H, especially in favoring relative expression of the maternal fish.
It is worth noting that a significant number of parental ELD pattern in DEG/DEM/DEP overlapped with the nonadditive pattern (383/1,122, 34% in DEG; 30/184, 16.3% in DEM; and 139/374, 37.2% in DEP). These observations suggest that the parental ELD pattern potentially contribute to heterosis in the hybrid, partly acting together with the nonadditive mechanism.
Interestingly, we found transgressive regulated proteins accounted for a larger proportion in the proteome DEP compared with those of DEG and DEM, especially in transgressive up-regulation (10 ϩ 102 ϩ 11), suggesting that the role of transgressive regulation is not to be ignored. Moreover, the significant KEGG pathways of NEP (5/12, 41.7%) were found to be overlapping for those of transgressive regulated proteins (Figs. 3 and 4). These situations imply that the ELD (dominance) of transcriptional and post-transcriptional levels changes to transgressive regulation (overdominance) of the protein level, which lays the molecular foundation for superior phenotype compared with two parents.
Comparative Analysis of KEGG Pathway Enrichment-To further explore the molecular mechanisms of H heterosis, KEGG pathway enrichment analyses of ELD-F/M genes/proteins were first conducted and then combined with those of nonadditivity and transgressive regulation for analysis. A total of 29, 16, 27, and 13 significant KEGG pathways (p Ͻ 0.05) were identified by ELD-F genes, ELD-M genes, ELD-F proteins, and ELD-M proteins, respectively (Figs. 4B-4G). By performing the KEGG pathway analyses together with ELD/ transgressive regulated genes/proteins, the biological process of "growth and development, feeding adjustment, metabolic capacity, and vitamin metabolism" tends to involve maternal fish, "digestive system" tends to involve maternal fish or dominant parents, and "cell proliferation" tends to involve paternal fish or the superior of the two parents. The "immune defense" and "antioxidant capacity" were frequently enriched by ELD/transgressive regulated genes/proteins ( Fig  4B-4G).
To further support ELD and the nonadditive mechanism for H heterosis, we found that these significant KEGG terms enriched from ELD-F/M genes/proteins were related to similar biological processes with NEG/NEP, such as immune defense, metabolism (vitamins, amino acids, lipids, and carbohydrates), digestion and absorption, and cell proliferation and development (Figs. 2, 4, and 5). The above results combined with dominance phenotypes of hybrid yellow catfish Huangyou-1 imply that both nonadditivity and parental ELD, especially high-parental ELD, play important roles in heterosis in H.

mRNA, miRNA and Protein Changes in Hybrid Catfish Heterosis
Correlation of DEM, DEG, and DEP Three Lines-Given that miRNAs negatively regulate the expression of their target mRNAs by target RNA cleavage, the expression patterns of miRNA target genes generally show an inverse correlation with those of miRNAs. Therefore, for most cases that involve target cleavage, the simple expectation is that when miRNAs are induced, their target mRNAs are reduced and vice versa. In consideration of the three comparison groups (M, F, and H), we conducted negative regulation analysis in 12 M-H-F expression patterns to dissect the possible miRNA-mRNA or miRNA-mRNA-protein interactions for heterosis, such as III(13)-IV(1,554)-IV(155), XϩXIϩXII(5 ϩ 7)-VIIϩVIIIϩIX(32 ϩ 170 ϩ 34)-VIIϩVIIIϩIX (10 ϩ 102 ϩ 11). The results showed that ELD-F "IV(46)-III(2,950)-III(281)" were the most common class in miRNA-mRNA interactions (Fig. 5A). There were 1612 negative miRNA-mRNA interactions, involving 46 DEM and 577 DEG in total, and 121 miRNA-mRNA-protein interactions, involving a total of 38 DEM and 31 DEG.
Data Extension and Validation of DEM, DEG and DEP by qRT-PCR and PRM-Although we believed that the data in this report had undergone rigorous statistical and bioinformatics analysis of three biological duplications to eliminate most of the possible errors, we nevertheless performed qRT-PCR and PRM analysis on 6 selected DEM/DEG/DEP to validate the mRNA-seq/miRNA-seq/TMT results (Fig. 5D). These miRNA/mRNA/proteins were manually selected as representatives for their potential roles in heterosis. These were found to be organism development (miR-27a-Camk2d1), metabolism (miR-30e-ACOT6 and miR-181a-Retsat), immune response (miR-221-Cfi and miR-130c-CHL1), and digestion and absorption (miR-140-Abcb11), according to the possible negative miRNA-mRNA-protein interactions on the basis of 12 male-hybrid-female (M-H-F) expression patterns.
The LC-PRM/MS analysis of PRM was performed on 7 peptide fragments of 6 target proteins from three groups of liver samples. All second-order mass spectrum matching figures of target protein candidate peptides were shown in supplemental Fig. S3. All original data of panorama were shown in supplemental Fig. S4. Data analyses of PRM quantitative Skyline in the target peptide section were shown in supplemental Table S12. The results indicate that quantitative data of target peptide fragments could be obtained from all nine samples (3 reps each of M, H, and F). Subsequently, relative quantitative analysis was performed on the target peptide fragments and the proteins through incorporation of heavy isotope-labeled peptide fragments. The results of qRT-PCR/ PRM revealed that most of these miRNA/mRNA/proteins share similar up/down-regulated expression tendencies with those from mRNA-seq/miRNA-seq/TMT data when compared with H (t test, p Ͻ 0.05). In addition, 60 genes involved in key biological processes of food intake regulation, growth and development, immune defense, antioxidant capacity, and flavor and nutrition were chosen for quantification by qRT-PCR to verify the transcriptome results and as additional evidence in the discussion (supplemental Fig. S5). Although there were some quantitative differences among these analytical platforms, the similarities between the RNA-seq/ miRNA-seq/TMT data and the qRT-PCR/PRM suggest that the RNA-seq data are reproducible and reliable. DISCUSSION The evolutionary force by hybridization results from the interactions between two different subgenomes in the newly formed hybrids. Thus, newly generated hybrids are good models to elucidate the heterosis and variation during hybridization. Our present results included the detailed information regarding parallel mRNA, miRNA, and protein expression levels in livers of P. fulvidraco, P. vachelli and a hybrid of both species, and we found a relatively large proportion of dominant regulatory changes. Given that miRNAs negatively regulate the expression of their target mRNAs by target RNA cleavage, providing the basis for protein changes to some extent, the expression patterns of miRNA target genes generally show an inverse correlation with those of miRNAs. Therefore, we predicted miRNA-mRNA-protein pairs and discuss these interactions and dominant regulatory changes here in the hope that the data will lay the foundation for the analysis of molecular regulation mechanism of heterosis.
In the present study, large amounts of NEG/NEM/NEP were found in the hybrid yellow catfish Huangyou-1 with several important physiological functions for heterosis (Fig. 2), which suggest unstable genomes and generation of the heterosis phenotype in progenies (2). Comparative genomics have indicated that during hybridization, the genome in progenies suffers a dynamic variation, including locus loss or duplication, which results in genome dominance, namely with genes from one subgenome having more loci and demonstrating greater expression (69,70). The genetic pattern for parental ELD and HEB genes has been found in fish allodiploid and allopolyploids, and these parental dominant models affecting the genotypes at the transcriptional level underlie heterosis in hybrids (12,13,15). Our result for H showed only 4.6% homoeologous genes (Table I, rows 13) with no bias of homoeolog expression, whereas the majority of homoeologous genes occurred with the HEB phenomenon and obtained either of the parental traits after the genome merger. Moreover, a large proportion of parental ELD pattern in H support the dominance hypothesis not only in transcriptional but also in post-transcriptional and protein levels (Fig. 4A). Meanwhile, the nonadditivity, HEB and ELD patterns were readily identified in different hybrid fishes, proving their widespread presence during hybridization (12,15). As can be seen from the color in the KEGG enrichment (Fig. 4), the number of highparental ELD pattern (types III and V) was significantly higher than that of the low-parental ELD pattern (types IV and VI) in transcriptome DEG and proteome DEP (p Ͻ 0.05, Fisher's exact test). We propose that the high parental expression of genes/proteins (growth, nutrition, feeding, and disease resistance) coupled with low parental miRNAs of the offspring are inherited from the mother or father, thus indicating that the offspring were enriched with the advantages of the father or mother and thereby obtained heterosis.
Based on the functionalities in the hybrid F1 of some plant species, the nonadditively expressed miRNAs identified were most likely key regulators that are critical for growth vigor and adaptation (18,71). Up to now, data on the small RNAs in the liver of fish is still limited, and information regarding small RNAs in fish heterosis is especially rare. Several nonadditive miRNAs were found in the hybrid cyprinidae (12). Of course, a few miRNAs were the same as the present study, such as miR-130, miR-27a, and miR-301, with well-characterized functions of immune responses and cell development. Some miR-NAs with functional features of metabolism, digestion, and absorption also contributed to H heterosis. Characterization and profiling of heterosis miRNAs in different fish is the beginning of a long road, and insights into miRNA-mRNA-protein regulatory networks (Fig. 5C) facilitate understanding of the fine-tuning of gene/protein expression at the post-transcriptional level (11).
In performing KEGG pathway analyses for NEG, NEP, and ELD/transgressive regulated genes/proteins, immune defense, metabolism (vitamins, amino acids, lipids, and carbohydrates), digestion and absorption, and cell proliferation and development were found to be the most frequently represented biological processes (Figs. 2,4,and 5). With regard to those heterosis phenotypes with faster growth and development, high feeding rate, and stronger immunity, we addressed our particular research question using pathway analysis to highlight the genes/proteins related to the functional clusters of (1) growth and development, (2) food intake regulation, (3) flavor and nutrition, and (4) immune defense and antioxidant capacity (Fig. 6). Most DEG, DEP, and their regulating miRNAs within these four functional clusters were validated using realtime PCR analysis and PRM (Figs. 5D and 6).
Growth and Development-In teleost fishes, as in other vertebrates, growth is regulated by insulin, the GH/IGF axis, local factors, and other endocrines (72,73). Insulin (INS) and Insulin-like growth factor (IGF) are major determinants of cell development and growth (74). In our results, the high parental ELD pattern genes (ELD-FϾM) of INS, Insulin receptor (INSR), and Insulin receptor substrate (IRS) suggested INS signal contributed to H rapid growth. IGF resembles INS in many respects, and they stem from a common precursor (75). Genes related to the GH/IGF axis such as IGF-1, IGFBP-1, IGFBP-2, IGFBP-4, and IGFBP-5a in the Hulong grouper liver were significantly up-regulated compared with one parent (4) and also increased IGF-1, IGFBP-1, and IGFBP-2 in hybrid cyprinidae liver (12). In this study, significantly abundant DEG of IGF-1 (ELD-M), IGFBP-3, IGFBP-2, and IGFBP-5 (ELD-FϾM) were also detected, which raised the GH/IGF axis, and H showed superior growth as compared with parents. These findings imply the significance of the GH/IGF axis which tends to be maternal fish for superior growth in hybrid fish.
The "Cancer" subclass was all enriched by ELD-M and transgressive regulated genes/proteins (Figs. 4D-4G). Although "Cancer" does not represent physiological and biochemical reactions for fish, the presence of more differentiated genes and proteins from the side was associated with cell proliferation and development (76). It is worth noting that "Rap1 signaling pathway," "PI3K-Akt signaling pathway," "HIF-1 signaling pathway," and "Cytokine-cytokine receptor interaction" may offer clues for transgressive regulated genes of cell proliferation and development involved in signal transduction, signaling molecules, and interaction. ELD-F genes associated with cell development were also found in hybrid cyprinidae (12).
In Hulong grouper brains, gonadotropin-releasing hormones were significantly up-regulated, inevitably affecting the downstream release of sex hormones (4,77). Previous studies have reported that androgen-related genes (e.g. 17␣-Methyltestosterone) could improve the GH/IGF axis and growth in P. fulvidraco (78,79). In H liver, estrogen-related genes were lower than the parents or tended to be linked to a parent with low expression, but androgen-related genes tended to be linked to a parent with high expression. Also observed was significant ELD-F KEGG enrichment of "Steroid biosynthesis" and "Steroid hormone biosynthesis" (also enriched by ELD-M proteins), which may be one of the causes of H rapid growth (Figs. 4B-4G).
Food Intake Regulation-In vertebrates, including fish, food intake regulation is the principal aspect of growth, which involves intricate networks of feeding adjustment, metabolic capacity, and the digestive system (80,81). For feeding adjustment, the up-regulated appetite-related genes adjusting the feeding rate of Hulong grouper have been reported (4). In our results, the presence of almost all genes/proteins reflecting an ELD-FϾM pattern among "Olfactory transduction" (DEG) and "Salivary secretion" (DEP) in H may have a stimulating effect on appetite contributing to food intake regulation.
Metabolic signals, arising from either the meal eaten or the changing pattern of substrate utilization, and hormone secretion in the inter-meal period may be monitored hepatically or centrally to contribute to the onset or termination of feeding (82). A strong metabolic capacity will increase an animal's food intake (83). In performing KEGG pathway analyses for ELD-F/M genes/proteins, fatty acid, carbohydrate, and amino acid metabolism were found to be the most frequently represented subclasses, which indicates that the vigorous anabolism and catabolism of lipids, saccharides, and proteins for H contributed to a strong metabolic capacity.
As for the digestive system, "Protein digestion and absorption" (all ELD-F gene/protein and transgressive regulated gene/protein), "Bile secretion" (ELD-M and transgressive regulated genes), "Gastric acid secretion" (ELD-F proteins), "Fat digestion and absorption," "Mineral absorption," "Pancreatic secretion" (transgressive regulated proteins), and "ABC transporters" (ELD-F) were significantly enriched, which is more evidence of the increased digestive ability for H. The above

mRNA, miRNA and Protein Changes in Hybrid Catfish Heterosis
findings imply the significance of food intake regulation for superior growth of the hybrid yellow catfish Huangyou-1.
Flavor and Nutrition-M, F, and H are edible commercial fish, and their flavor and nutrition are important aspects of their commercial value. Vitamins are necessary for maintaining health. For instance, thiamine (Vitamin B1) not only is an important substance affecting the flavor and taste of meat but also promotes gastrointestinal motility and improves human appetite (84). In this study, the significantly enriched KEGG by the up-regulated genes for "Thiamine metabolism" (ELD-F), "Retinol metabolism" (Vitamin A; ELD-F; also enriched by ELD-M proteins), "Pantothenate and CoA biosynthesis" (ELD-M), and "One carbon pool by folate" (transgressive regulation) as well as significant KEGG by ELD-FϾM proteins of "Riboflavin (B2) metabolism," "Vitamin B6 metabolism," and "One carbon pool by folate" suggest that the metabolic capacity of vitamins A, B1, B2, B6, pantothenic acid (B5), and folic acid (B9) in H liver were higher than in the parents or tend to come from a parent with high expression.
Essential fatty acids (e.g. linoleic acid) are as important as vitamins and proteins, and human clinical trials have proven polyunsaturated fatty acids can relieve acute and chronic inflammatory diseases. The significant KEGG enrichment of "Biosynthesis of unsaturated fatty acids" (ELD-M) and "Linoleic acid metabolism" (ELD-F) with high parental ELD pattern genes implies that the biosynthesis of unsaturated fatty acids such as linoleic acid are linked to a parent with high expression. In addition, "Fatty acid biosynthesis" (ELD-F proteins) and "Taurine and hypotaurine metabolism" (ELD-F and transgressive regulated genes) were significantly enriched, which suggests stronger ability to metabolize nutrients.
Immune Defense and Antioxidant Capacity-Hybrid fish are generally more resistant to disease and oxidation (85), allowing them to resist bacteria, exogenous substances, and other factors, thereby increasing their survival rates (86). KEGG pathway analyses for ELD and transgressive regulated genes/ proteins showed "Human diseases," "Immune system," and "Xenobiotics biodegradation and metabolism" to be the three most frequently represented subclasses for immune resistance. Although "Human diseases" does not represent physiological and biochemical reactions for fish, the presence of more differentiated genes and proteins from the side has been associated with resistance to disease (87).
In the immune systems, cell adhesion molecules are glycoprotein expressed on the cell surface and play a critical role in a wide array of biologic processes, including immune response and inflammation. The main concentrations of proteins associated with drug transport were ABCA, ABCB, ABCC, and ABCG (88). It is noteworthy that "Cell adhesion molecules (CAMs)" and "Phagosome" were all enriched by ELD-F/M genes and transgressive regulated proteins, and the differential genes are mainly related to ABCA, ABCB, ABCC, showing that H liver has a better capacity to recognize specific receptors and phagocytic ability. The complementary system is a proteolytic cascade in blood plasma and a mediator of innate immunity, a nonspecific defense mechanism against pathogens (89). The enrichment of "Complement and coagulation cascades" (all ELD-F/M genes/proteins) exerted the high parental ELD pattern contributed to a stronger innate immunity, which is the same with the expression pattern of complement activity for C. carpio ϫ C. gibelio (90). We found that the most up-regulated expressed genes/proteins of the immune pathway, just like "Hematopoietic cell lineage" (ELD-M genes/proteins), "B cell receptor signaling pathway" (ELD-M genes), and "Chemokine signaling pathway" (transgressive regulated genes), more favored either the paternal parent or both the parents, apart from the common KEGG enrichment pathways, which indicated that H has a better immune system than the parent.
The "Xenobiotics biodegradation and metabolism" subclass highlights the detoxification of metabolic by-products and oxidative stress in the presence of more oxidation resistance in H, as "Drug metabolism" and "Metabolism of xenobiotics by cytochrome P450" were frequently enriched by ELD genes/proteins, as well as resistance to disease and stress resistance, similar to the hybrid cyprinidae (12). Excessive reactive oxygen species (ROS) can lead to oxidative stress, and ROS promptly attacks cellular components. The effective and rapid elimination of ROS is accomplished by antioxidant defense systems such as superoxide dismutase (SOD), catalase (CAT), and glutathione peroxidase (GPx), which can result in the detoxification of ROS through simple redox reactions. The function of these enzymes is strongly associated with increased levels of glutathione S-transferase (GST) and glutathione (91). In this study, the significantly abundant genes GPx (ELD-M), SOD3 (ELD-MϾF), CAT, GSTA4, GSTT1, GSTCD, GSTP1, GSTO1, GGCT (Gamma-glutamylcyclotransferase), and ANPEP (Aminopeptidase) (ELD-FϾM) and ELD proteins of GPx, glutathione synthetase (GSS), GSTA1, GSTO1, GSTT1, and GSTT3 were detected, which raised oxidation-reduction. "Glutathione metabolism" (ELD-F genes/proteins) and "Peroxisome" (ELD-M proteins) were also observably enriched. This demonstrates that H combined the qualities of parents with higher antioxidant capacity than the parents.
Conclusions and Perspective-We provide a good case study for comprehensive analysis of mRNA/miRNA/protein expression and profiling of non-model species. Our data show that the nonadditivity, HEB, and ELD patterns were readily identified in the hybrid yellow catfish Huangyou-1, providing evidence of the widespread presence of dominant models during hybridization. Furthermore, some of the diverse key pathways studied, including immune defense, metabolism, digestion and absorption, and cell proliferation and development, helped us articulate the vital mechanisms involved in the generation of heterosis phenotype in progenies (Fig. 6). Our data suggest a potential miRNA-moderated regulatory mechanism for high parental gene expression as well as protein expression as the basis for heterosis phenotypes in hybridization. The above results provide new and important information about the molecular mechanisms of heterosis, representing a significant step toward a more complete elucidation of this phenomenon.
Despite the apparent conceptual simplicity of the construction of miRNA-mRNA-protein regulatory networks, the genetic randomness of hybrids, gene expression noise, interaction of three levels, and so on, are not yet well-characterized. Thus, the combination of multilevel high throughput deep sequencing datasets (e.g. chromatin immunoprecipitation, methylation or degradome sequencing data) with bioinformatics analysis could serve as a powerful tool for building better networkbased molecular models for predicting, testing, and identifying robust miRNA-mRNA-protein pairs and the main biological processes for heterosis. Taken as a whole, although this study cannot completely account for all molecular mechanisms in fish heterosis, we hope that it provides basic data in terms of analytical methodology and new biological and mechanistic insight for future research.

DATA AVAILABILITY
The datasets generated and analyzed during the current study are available in repositories as follows: the high-throughput sequencing data accession codes: GSE118627 (mRNA-seq and miRNA-seq dataset) (https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?accϭGSE118627). TMT proteomics data and PRM data related to this study has been publicly available on iProX (www.iprox.org) with ID IPX0001286000 and IPX0001362000.