Integrated Analysis of mRNA and miRNA Changes in Two Haliotis diversicolor Genotypes and Their Hybrid

Heterosis is a widely distributed phenomenon in mollusks. It is vital in aquaculture by bringing beneficial traits into hybrids. People have utilized the heterosis theory in aquaculture for years. However, the molecular basis of heterosis remains elusive. Evident growth and survival heterosis were shown in the hybrid (“Dongyou-1”) of two Haliotis diversicolor geographic genotypes (Japan and Taiwan). To explore the molecular basis underlying the hybrid abalone’s heterosis, we conducted comparative mRNA and miRNA transcriptional analysis in the hybrid and parental genotypes. Differentially expression analysis identified 5,562 differentially expressed genes (DEGs) and 102 differentially expressed miRNAs (DEMs) between the three genotypes. 1,789 DEGs and 71 DEMs were found to be non-additively expressed in the hybrid. Meanwhile, both the expression level dominance pattern (ELD) and expression level overdominance pattern (ELOD) were found in the DEGs and DEMs, showing the existence of dominance and overdominance models in the hybrid’s transcriptome and post-transcriptional regulation. Functional analysis showed the non-additively expressed genes, ELD genes, and ELOD genes were significantly enriched in growth, immunity, and stress response related pathways, while some of the pathways were regulated by the mRNA-miRNA interactions. The expression levels of FGF, C1Q, HC, CAT, SEGPX, and MGST were significantly up-regulated in the hybrid compared to the middle parent value. In conclusion, we identified the existence of non-additivity, dominance, and overdominance models in the transcriptome and miRNAome of the H. diversicolor hybrid; these models facilitate the advantageous parental alleles’ integration into the hybrid, contributing to the hybrid’s growth and survival heterosis.


INTRODUCTION
Heterosis, or hybrid vigor, refers to the phenomenon that hybrids often show advantages in growth, fitness, or resistance relative to the parents (Chen, 2013). Heterosis is widespread in plants and animals, such as Arabidopsis, maize, rice, cotton, mice, donkey, chicken, dog, oyster, abalone, and mussel (Chen, 2013;McKeown et al., 2013). Farmers and researchers have utilized the heterosis theory to improve agricultural species' genotypes and phenotypes for years (Goff, 2011;Veitia and Vaiman, 2011;Baranwal et al., 2012;McKeown et al., 2013). Heterosis is also common in the interspecific hybrids of abalones. For example, H. discus hannai ♀ × H. discus ♂ and H. discus hannai ♀ × H. fulgens were superior in growth or viability than the parents (Inoue et al., 1986;You et al., 2015). The reciprocal hybrids of H. gigantea and H. discus hannai were more resistant to heat and bacteria than the parents (Luo et al., 2010;Liang et al., 2014Liang et al., , 2018. Besides, heterosis also exists in the intraspecific hybrids of abalones. For instance, the hybrid of the H. discus hannai China genotype and Japan genotype showed more considerable shell length, shell width, wet weight, and a higher survival rate (17.9, 22.1, 61.9, and 80%, respectively) compared to the China genotype (Zhang et al., 2004). Some abalone hybrids have been successfully applied in abalone aquaculture due to their superior traits Cook, 2019).
Although abalone heterosis is important in aquaculture, its molecular mechanisms remain unclear. Traditional analyses of molluscan species' heterosis mainly focused on the relationship between isozyme heterozygosity and heterosis, or heterosisrelated genes and proteins. For example, a positive correlation between organisms' isozyme heterozygosity and growth rate was identified in mussels and oysters (Gentili and Beaumont, 1988;Hedgecock et al., 1995Hedgecock et al., , 1996. The hybrid H. discus hannai ♀ × H. gigantea ♂ had higher activities of immunityrelated enzymes (superoxide dismutase, acid phosphatase, alkaline phosphatase, and myeloperoxidase) and higher mRNA level of heat shock protein 70 gene upon heat stress compared to the parents (Liang et al., 2014). Similarly, the hybrid Haliotis discus hannai ♀ × H. fulgens ♂ showed higher phenoloxidase activities, superoxide dismutase, and higher mRNA levels of several immunity-related genes compared to the maternal genotype (Shen et al., 2020). These studies give us a deeper understanding of molluscan heterosis by identifying heterosis-related molecules. However, to have a comprehensive understanding of molluscan heterosis, we need to conduct further genome-scale studies of heterosis using omic tools.
Transcriptome sequencing is the most frequently used omic tool in heterosis studies. It has been utilized in the heterosis studies of many species such as Arabidopsis, rice, wheat, corn, fish, and oyster (Hedgecock et al., 2007;Goff, 2011;Chen, 2013;Gao et al., 2013). Although different species have different parenthybrid gene expression variations, some common trends could still be concluded. The parental alleles' expression in the hybrid could be divided into several patterns: additivity, non-additivity, dominance, and overdominance (Baranwal et al., 2012). For the non-additivity, dominance, and overdominance patterns, the heterozygote's gene expression value is unequal to the middle parental expression value (MPV) (Chen, 2010); in other words, a new allelic expression pattern is generated in the hybrid. Hence, the non-additivity, dominance, and overdominance patterns are vital to the heterosis generation (Springer and Stupar, 2007;Fujimoto et al., 2012). The differentially expressed genes (DEGs) among parents and hybrids are functionally enriched in various biological pathways such as "cell growing and aging", "stress response", and "signal transduction" (Chen, 2013). These pathways are closely related to organisms' growth, fitness, and resistance. Besides, the role of microRNAs (miRNAs) in the heterosis generation obtains much attention in recent years. miRNAs can negatively regulate gene expression through epistasis (Sha et al., 2014). The miRNAs expression variation between hybrids and parents has been investigated in many species such as Arabidopsis, wheat, and corn (Greaves et al., 2012;Chen, 2013). A common finding is the number of non-additively expressed miRNAs (NEMs) in hybrids positively correlates with the parents' genetic distance. The NEMs could further cause allele expression changes in hybrids (Groszmann et al., 2011(Groszmann et al., , 2013Ng et al., 2012). To conclude, the transcriptional analysis could help us understand how parental mRNAs and miRNAs are integrated and expressed in hybrids. Meanwhile, the joint analysis of transcriptional data and physiological data could help us identify essential genes or pathways related to the heterosis generation.
Haliotis diversicolor is an important economic abalone species cultured in China's continental coasts. This species was initially introduced from Taiwan and then successively cultured for several generations. However, due to frequent inbreeding, the TW's survival rate continuously declined. To solve this problem, we introduced several different geographical genotypes of H. diversicolor to hybridize with TW. The hybrid genotypes showed significantly different growth and survival traits compared to the parents. Within the hybrid genotypes, the hybrid of H. diversicolor Taiwan genotype and Japan genotype (JP) showed the most desirable traits; it was TWlike in growth and JP-like in survival . This hybrid (named DY) has been widely used in China's abalone aquaculture. However, we know little about the molecular mechanisms underlying DY's heterosis. In this study, we conducted comparative transcriptional analysis in DY and its parents (JP and TW) to compare the mRNAs and miRNAs differences between the three genotypes. We wish to identify the important biological molecules and functions determining DY's heterosis and investigate the miRNA-mRNA interactions in the abalone's heterosis generation.

Abalones and Phenotype Measurements
Abalones of the three H. diversicolor genotypes (JP, TW, and DY) were used. The abalone breeding and cultivation were performed at "Fuda" abalone company in Jinjiang City, Fujian Province, China. The artificial fertilization of JP ♀ × JP ♂, TW ♀ × JP ♂, and TW ♀ × TW ♂ was conducted in October 2015. The obtained abalone seeds were labeled and mixed cultured under the same conditions (reared in running seawater; fed with fresh alga Gracilaria confervoides; salinity 33 ; pH 8.0) for 2 years. The three genotypes' survival rates were recorded every month from October 2016 to March 2017. The shell length and wet weight of 100 abalones from each genotype were measured in October 2017.

Acclimation
For each genotype, 30 abalones with a normal appearance and good vitality were used for the following experiments. The abalones were acclimated under the controlled laboratory conditions (reared in still seawater with aeration, the seawater was replaced every day; fed with fresh alga Gracilaria confervoides; temperature 24 ± 1 • C; salinity 33 ; pH 8.0) for 14 days. No abalone was dead during the acclimation.

RNA Extraction
For each genotype, ten abalones were used for the transcriptome and miRNAome sequencing. Total RNA was extracted from each abalone's muscle using the TRIpure reagents (Invitrogen, United States). The RNA quality was tested using 1% agarose gels and Agilent 2,100 Bioanalyzer (Agilent Technologies, United States). The ten RNA samples from each genotype were equally mixed together (400 ng RNA from each sample) and then equally divided into two RNA pools (2 µg RNA for each pool). The two RNA pools were then used to construct the cDNA library and small RNA library separately.

cDNA Library Construction and mRNAs Sequencing
For each genotype, one RNA pool, as described below, was used for the poly (A) + mRNA selection using oligo (dT) magnetic beads (Invitrogen, United States). The cDNA library was constructed using the ScriptSeq kit (Epicentre, United States). For the mRNAs sequencing, a total of three cDNA libraries were constructed. The libraries were amplified by 12-15 cycles of PCR and then sequenced in two lanes on the HiSeq 2,500 sequencer (Illumina, United States) at BGI-tech company (Shenzhen, China).

De novo Assembly and Functional Annotation
The generated raw reads were first filtered by removing the adaptors, ambiguous reads (those with ≥10% ambiguous nucleotides), and low-quality reads (those containing ≥20% bases whose quality scores ≤10) then assembled using Trinity (Grabherr et al., 2011), with min_kmer_cov set to 2 and other parameters set to defaults. The obtained contigs and unigenes were then aligned to the non-redundant (NR) protein database, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) using KAAS (Mao et al., 2005;Moriya et al., 2007), with an E-value cut-off of 10 −5 .

Differentially Expression Analysis
The read counts were normalized as FPKM (Fragments Per Kilobase of exon model per Million mapped reads) values. Then statistical analyses were applied on FPKM of genes in the parents and hybrid using DEseq 1 . DEGs between any two populations were defined to have a significant expression level against each other (fold-change ≥1.5, p-value ≤ 0.01). Nonadditively expressed genes (NEG) and NEMs in the hybrid were 1 http://www.bioconductor.org/packages/release/bioc/html/DESeq/ defined to have a significantly different expression level against MPV (fold-change ≥1.3, p-value ≤ 0.01).

Small RNAs Sequencing
As described before, the three RNA pools (one from each genotype) were used to construct three small RNAs libraries for small RNAs sequencing according to the manufactures of TruSeq Small RNA Sample Prep Kits (Illumina, United States). For each library, 1 µg of total RNA was used. The three libraries were then sequenced on the HiSeq 2500 sequencer (Illumina, United States) at BGI-tech company (Shenzhen, China).

Bioinformatic Analysis of Small RNAs Sequencing Data
The 49 nt sequence tags from Hiseq sequencing went through the data cleaning analysis to get credible clean tags without adaptor and unknown nucleotides (>10%). Only 18-32 nt sequences were used for further study. Then the length distribution of the clean tags and common and specific sequences between samples were summarized. Then clean tags were aligned to the Rfam database 2 using bowtie (Langmead et al., 2009) to identify miRNA sequences. The known miRNA sequences were obtained by aligning the clean tags to the miRNase database 3 by BLAST 4 , and the novel miRNA sequences were predicted using miRDeep2 (Friedländer et al., 2012).

Differentially Expression Analysis and Target Genes Prediction
Differentially expressed miRNAs (DEMs) between the parents and hybrid were analyzed as the method described in "2.7 Differentially expression analysis". The target genes of the DEMs were identified by targetscan (García et al., 2011), miRanda 5 , PITA 6 , and RNAhybrid (Rehmsmeier et al., 2004). The overlapping target genes predicted by the software were kept for further analysis.

Twelve Expression Patterns of DEGs and DEMs
To further dissect parental alleles' expression in DY, the DEGs and DEMs between the parents and hybrid were binned into 12 categories according to the previous studies (Rapp et al., 2009;Zhang et al., 2019). To investigate the possible miRNA-mRNA interactions for abalone heterosis, we achieved the negative interaction between the DEMs and DEGs by integrating DEGs' expression profiles and DEMs with the DEMs targeting genes information as described before .

Functional Analysis of DEGs and DEMs
Gene Ontology analysis of the DEGs and DEMs (target genes of the DEMs) was conducted by Blast2GO 7 with a significant threshold of FDR < 0.05. KEGG pathway analysis of the DEGs and DEMs was performed through the hypergeometric test using the corrected p-value < 0.05 by KOBAS (Mao et al., 2005).

Expression Profiles of Growth, Immune, and Stress Response Genes
Based on NR annotation data and previous studies Wang et al., 2018b;Bouallegui, 2019), 30 growth-related genes, 52 immunity-related genes, and 29 stress response genes were identified in the three genotypes. For comparison, each genotype's gene expression level was normalized by comparing to the MPV. Besides, each gene group (genes with the same functions) was clustered between the three genotypes based on the Pearson correlation algorithm.

Data Visualization
Heat maps of gene expression were drawn using Morpheus 8 while the genes were clustered based on the Pearson correlation algorithm. Venn diagrams were drawn by Venny 9 .

Real-Time PCR Validation
To validate the RNA-seq results, we assessed the expression levels of six DEGs (unigene 34268, unigene 47208, unigene 47561, unigene 32067, unigene 19086, and unigene 49863) and six DEMs (miR-4755-5p, miR-6317, miR-7255-5p, miR-4326, miR-2836, and miR-263b) by real-time PCR, according to the method described below (Liang et al., 2018). Briefly, the realtime PCR experiment was conducted in a 7500 fast qPCR system (ABI, United States) using specific qPCR reagents (Thermo, United States). The cycling parameters used were: 95 • C for 7 min, 35 cycles at 95 • C for 20 s, and 60 • C for 1 min. The relative mRNA level of target genes was calculated based on the Ct values of this gene and β-actin normalized to that of the cDNA standard.

Phenotypes of the Three H. diversicolor Genotypes
The three genotypes' survival rates from October 2016 to March 2017 are shown in Supplementary Figure 1. During the 6 months, the JP or DY's survival rate was significantly higher than that of TW. By March 2017, JP, DY, and TW's survival rates were 93.56, 95.56, and 37.56%, separately. The DY's shell length and wet weight were similar to that of TW, higher than that of JP (Supplementary Table 1). These results suggest DY has JP-like resistance and TW-like growth ability.

Summary of Transcriptomes and miRNAomes
After the sequencing quality control, an average of 54 million raw reads and 52 million clean reads were obtained for each transcriptome sample. The average Q20 percentage of the samples was 97.01% (Supplementary Table 2). After assembling,85,979,75,220,and 72,936 unigenes were generated (with an average length of 750, 667, and 672 nt) in JP, DY, and TW, respectively. The N50 of all unigenes was 2,166, and the mean unigene length was 1,075 nt (Supplementary Table 3).
An average of 11,455 thousand raw reads and 11,318 thousand clean reads were obtained for each miRNAome sample. Most of the reads were 20-30 bp in length (Supplementary Table 4). After alignment and target genes prediction, 3,186, 2,903, and 3,142 known miRNAs were identified in JP, DY, and TW, corresponding to 68,633, 68,623, and 68,632 target genes, respectively. 20, 14, and 9 novel miRNAs were identified in JP, DY and TW, corresponding to 34,363, 31,985, and 26,692 target genes, respectively (Supplementary Table 5).

DEGs and DEMs
Differentially expression analysis identified 3,661, 3,048, 1,778 DEGs and 65, 63, and 68 DEMs between JP and TW, DY and JP, and DY and TW, separately (Figure 1). The DEGs numbers suggest a distant gene expression pattern between the two parental genotypes, while DY's gene expression pattern was more like TW other than JP. The number of DEMs was almost the same between any two genotypes. Besides, there were 99 overlapping DEGs and one overlapping DEM between the three genotypes.

Functional Analysis of NEGs
The NEGs were mainly enriched in "biological process" category GO terms. The "development process", "chromosome", and "phosphoric ester hydrolase activity" were the most enriched GO terms in the biological process, cellular component, and molecular function categories, respectively (Supplementary Figure 2). KEGG pathway analysis of the NEGs identified 35 significantly enriched pathways (p < 0.05) (Figure 3). Most of the enriched pathways were related to "organismal systems (digestive, endocrine, excretory, immune, sensory, and circulatory systems)", "diseases (immune, infectious, neurodegenerative, and cardiovascular diseases)", and "metabolisms (amino acid, xenobiotics, and carbohydrate metabolisms)". The enrichment results suggest the hybrid is different from the parents' average performance (MPV) in growth, immunity, and stress-response-related biological functions; these differences may contribute to the DY's superior phenotypic traits.

Twelve Expression Patterns of DEGs and Functional Analysis
The 12 categories are shown in Figure 4: (1) Additivity (I and II); (2) Expression level dominance by JP (ELD-J, III and IV); (3) Expression level dominance by TW (ELD-T, V and VI); (4) Expression level underdominance (ELUD, VII, VIII, and IX); and (5) Expression level overdominance (ELOD, X, XI, and XII). Notably, a large number of ELD patterns were identified, with the proportion of 52.28% (1,226/2,345) in DEGs and 61.64% (45/73) in DEMs. The second large category was the ELOD pattern, with the proportion of 8.70% (204/2,345) in DEGs and 24.66% (18/73) in DEMs. Besides, the number of ELD-T patterns was significantly higher than that of ELD-J patterns in DEGs, while in DEMs the two patterns had similar numbers. These findings indicate the widespread of ELD and ELOD patterns in the hybrid's transcriptome and miRNAome.
Gene Ontology analysis showed the ELD-J genes were mostly enriched in "biological adhesion" and "gland development" (biological processes); "extracellular region" and "late endosome" (cellular components); "isomerase activity" and "binding-related activities" (molecular functions). The ELD-T genes were most enriched in "development process" and "organ development" (biological processes), "non-membrane-bounded organelle" (cellular components), "binding-related activities" (molecular functions). The ELOD genes were mainly enriched in "transportrelated processes" (biological processes), "membrane-related components" (cellular components), and "transmembranerelated functions" (molecular functions) ( Supplementary  Figures 3-6). These results suggest the ELD-J and ELOD genes are more related to cellular interaction and signaling transduction functions, while the ELD-T genes are more related to growth and development functions.
KEGG pathways analysis identified 11, 11, 7, and 18 significant (p < 0.05) pathways from the ELD-J, ELD-T, ELOD, and ELUD genes separately (Figure 5 and Supplementary  Figure 7). Pathways enriched from the ELD-J genes were mainly related to "diseases (immune disease)", "organismal systems (digestive, nervous, and circulatory systems)", and "environment information processes". Pathways enriched from ELD-T genes were mostly related to "diseases (cardiovascular, infectious, and neurodegenerative diseases)" and "metabolism (carbohydrate and lipid metabolisms)". Pathways enriched from the ELOD genes were mainly related to "metabolisms (amino acid, cofactors, and vitamins metabolisms)", "organismal systems (digestive, endocrine, immune, and excretory systems)", and "environmental information processing (signal transduction, signaling molecules, and interaction)". The above pathways' functions are closely related to organisms' growth, immune response, and organism's reaction to external environments. Hence, the ELD and ELOD genes are crucial to the DY's heterosis generation.
Interestingly, by comparing the pathways enriched from the ELD and ELOD genes with those enriched from the NEGs, we found 63.67% (7/11) of the ELD-J pathways, 45.46% (5/11) of the ELD-T pathways, and 50.00% (9/18) of the ELOD pathways were overlapped with the NEGs pathways (Supplementary Datasheet 1). The overlapped ELD pathways were mostly related to "immune disease", "nervous and circulatory systems", and "signaling molecules and interaction", while the overlapped ELOD pathways were mainly related to "amino acid metabolism" and "organismal systems (digestive, endocrine, immune, and excretory systems)". The above results suggest the ELD and ELOD patterns partially act together with the non-additivity pattern in altering the hybrid's growth, immunity, and stress response related functions, which might be a vital mechanism for the emergency of abalone heterosis.

Interactions Between DEGs and DEMs
As Figure 6 shows, the most common class in the miRNA-mRNA interactions was ELD-J "IV(19)-III(173), " involving 19 DEMs and 115 DEGs. Besides, 50.32% of the ELD-J genes Gene Ontology analysis showed the ELD-J pairs genes were mostly enriched in "cell adhesion" and "biological adhesion" (biological processes); "nucleolus" (cellular components); "DNA binding" and "transcription regulator activity" (molecular functions). The ELD-T pairs genes were most enriched in FIGURE 3 | KEGG pathway enrichment of the NEGs. This figure consists of the pathways, the second-class hierarchy, and the first-class hierarchy. The first-class hierarchy and second-class hierarchy rank from top to bottom according to the number of pathways included in the hierarchy. The pathways rank from top to bottom according to the enrichment significance of each pathway. Each bar represents the enrichment count of transcripts in a KEGG pathway. *EIP, environmental information processing; CP, cellular processes; GIP, genetic information processing. "development process" and "organ development" (biological processes), "organelle" (cellular components), "bindingrelated activities" (molecular functions). The ELOD/ELUD pairs genes were mainly enriched in "DNA integration" (biological processes), "respiratory-related complexes" (cellular components), and "oxidoreductase activity (NADPH)" (molecular functions) (Supplementary Figures 8-10). The enriched GO terms of ELD, ELOD, and ELUD pairs genes are quite similar to that of the ELD, ELOD, and ELUD genes.

Expression of Growth, Immunity, and Stress Response Related Genes
As is shown in Figures 8-10, most of the growth-related genes showed higher expression levels in DY and TW other than JP. Notably, the fibroblast growth factor receptor gene (FGFR) was overdominant expressed in DY. These results confirm the existence of dominance and overdominance models in DY's growth-related genes expression.
DY had a significant maternal expression of most immunityrelated genes (13 out of the 20 gene groups) and a significant paternal expression of most stress response genes (four out of the seven gene groups). This indicates the similar survival phenotype of JP and DY (Supplementary Figure 1) might be due to their similar expression of stress response molecules, other than immunity-related molecules. Besides, we also found DY has significantly higher expression levels of complement C1 gene (C1QI2), hemocyanin gene (HC-2), catalase gene (CAT), glutathione peroxidase gene (SEGPX), and glutathione s-transferase gene (MGST1) than MPV. These genes might also be the important molecules related to DY's heterosis generation.

DISCUSSION
In recent years, transcriptional analysis has been conducted in various aquatic species to compare the allelic expression variations between the parents and hybrids, such as brook charr (Bougas et al., 2010(Bougas et al., , 2013, putterfish (Gao et al., 2013), grouper (Sun et al., 2016a,b), crucian (Ren et al., 2016), and bream (Zheng et al., 2019). However, the transcriptional analysis of molluscan heterosis was only reported in the oysters (Hedgecock et al., 2007;Yang et al., 2018) and the sea cucumber (Wang et al., 2018a). In this study, we compared the growth and survival phenotypes of two H. diversicolor genotypes and the hybrid, analyzed the transcriptomes and miRNAomes of the three genotypes, and predicted the miRNA-mRNA interactions. We wish to explore the molecular mechanisms of the H. diversicolor heterosis through this study and lay the foundation for molluscan heterosis analysis.
DY integrated profitable maternal growth and paternal advantageous survival phenotypes (Supplementary Figure 1 and Supplementary Table 1). In our previous study, DY was JP-like in low-temperature and air exposure resistance (You et al., 2019). The hybrids' inheritance of parental advantages is not uncommon in mollusks. For example, the F 1 hybrid of H. discus hannai and H. kamtschatkana was fast in growth like H. discus hannai and strong in resistance like H. kamtschatkana (Hoshikawa et al., 1998). The F 1 hybrid of H. discus hannai and H. gigantea had maternal growth superiority and paternal resistance superiority (Luo et al., 2010;Liang et al., 2014Liang et al., , 2018. A similar phenomenon was also observed in the oyster (Guo et al., 2008), mussle (Beaumont et al., 2004), and scallop (Zheng et al., 2006). The hybrids' vigorous phenotypes may result from the interactions between two distant parental subgenomes.
The ELD and ELOD expression patterns have been identified in the hybrids of oyster (Hedgecock et al., 2007), pearl oyster (Yang et al., 2018), and various fish species (Bougas et al., 2013;Gao et al., 2013;Ren et al., 2016;Sun et al., 2016a;Zhang et al., 2019). In this study, large amounts of ELD and ELOD genes were also identified in the hybrid abalone (Figure 4). Interestingly, the ELD-J genes were mostly high-parent patterned (type III) rather than low-parent patterned (type IV), while the ELD-T genes were mostly low-parent patterned (type VI). In the yellow catfish and the Cyprinidae's hybrids, both the maternal and paternal ELD genes were high-parent patterned (Zhou et al., 2015;Zhang et al., 2019). The authors of the two studies propose the hybrids may inherit advantageous molecules from both parents to obtain heterosis. In our case, the hybrid abalone may inherit more advantageous alleles from the paternal genotype other than the maternal genotype. Besides, the number of ELOD genes was much higher than ELUD genes, indicating the overdominance models are more common in the hybrid than the underdominance models. This phenomenon agrees with the yellow catfish and Cyprinidae studies. The overdominance model is an essential mechanism of transgressive heterosis by creating a more robust heterozygous gene expression pattern (Baranwal et al., 2012;Chen, 2013). Moreover, the ELD and ELOD patterns were also found in the DEMs (Figure 6), indicating the dominance and overdominance models exist in the transcriptome and post-transcriptional regulation of the hybrid abalone.
Most of the ELD genes, ELOD genes, ELD pairs genes, and ELOD pairs genes were functionally related to metabolisms, organismal developments, diseases, and environmental information processing (Figure 5). Considering the DY's heterosis in growth and survival, we could discuss these genes by three functional categories: (1) Growth and development, (2) Immune defense, and (3) Stress response.
(1) Growth and development The growth-related pathways were mostly enriched from ELD-T and ELOD genes, including metabolism, growth, digestion related pathways. ELD-J genes were also enriched in two metabolism pathways even though the enrichment significances were relatively small (Figure 5). These pathways are crucial for organisms to keep normal growth functions (Gao et al., 2013). In other aquatic species, ELD or ELOD genes' functions are also frequently related to growth Wang et al., 2018b). Hence, hybrids' growth heterosis could result from their similar growth-related functions with the FIGURE 10 | Relative expression of stress response genes in the three H. diversicolor genotypes. * represents the significant difference of the gene expression level compared to the MPV. more robust progeny, or their transgressive functions beyond both parents. In this study, DY was TW-dominant in the carbohydrate and lipid metabolisms; overdominant in the protein and lipid digestions, amino acids, and vitamins metabolisms. These dominant or overdominant biological functions may work together in forming DY's growth phenotypes. Besides, the "carbohydrate metabolism, " "cell growth and death, " and "amino acid metabolism" pathways were also enriched from the ELD-T pairs and ELOD/ELUD pairs genes (Figure 7), suggesting the existence of miRNA regulations in these biological functions. The overlapping of ELD gene pathways and ELD pairs genes pathways was also observed in the hybrid yellow catfish .
The 30 growth-related genes were mostly ELD-T, or ELOD expressed in DY (Figure 8). The dominant or overdominant expression of growth-related genes was identified in almost all the comparative transcriptional analyzes between hybrids and parents (Di et al., 2013(Di et al., , 2015Sun et al., 2016a;Zhang et al., 2017;Wang et al., 2018a;Zheng et al., 2019). Functions of the ELD and ELOD genes differed within different species, showing the different heterosis generation mechanisms. Notably, the FGFR gene was overdominant expressed in DY. The FGFR family genes are tyrosine kinase receptors controlling the functions of fibroblast growth factor (FGF) genes by binding and activating processes (Dailey et al., 2005;Beenken and Mohammadi, 2009). Four FGFR transcripts were found to be highly expressed in DY, indicating DY's superior FGF-related functions. The proteomic analysis of the three H. diversicolor genotypes also found several growth-related proteins to be overdominant expressed (ATP synthase subunit, fructose-bisphosphate aldolase, and triosephosphate isomerase) or partial dominant expressed (enolase, arginine kinase, and taurine dehydrogenase) in the hybrid (Di et al., 2013). These results suggest the ELD and ELOD expression patterns of growth-related molecules exist in the hybrid's transcriptional and proteomic levels.
(2) Immune defense In molluscan heterosis studies, the immunity-related molecules draw less attention than the growth-related molecules. However, the immune ability is also crucial to hybrids by improving the organisms' resistance (Miller et al., 2015;Ma et al., 2018;Zhao et al., 2019). In our study, the diseases related pathways were mainly enriched from the ELD-J and ELD-T genes, suggesting the integration of parental immune characteristics into the hybrid. The ELOD genes were also enriched in several immunity-related pathways, indicating the DY's transgressive performance in some immune functions. However, DY's survival phenotype highly conformed with JP and differed from TW (Supplementary Figure 1). Hence, DY's immunity-related molecules' expression may not be enough to explain its survival heterosis toward TW. Besides, the "infectious disease" and "immune disease" pathways were also enriched from the ELD pairs and ELOD/ELUD pairs genes, suggesting the role of miRNA regulations in some immunity-related pathways during hybridization.
The mollusks' immunity-related molecules could be divided into three functional categories: pattern recognition receptors, signaling molecules, and effector molecules (Wang et al., 2013a;Song et al., 2015;Bouallegui, 2019). Of the 52 immunityrelated genes investigated, DY and TW had lower expression levels of most genes (11 out of the 19 gene families) than JP. However, DY showed an overdominant expression pattern in two genes, C1QI2 and HC-2 (Figure 9). The C1Q gene is a classic peptidoglycan recognition protein gene containing the versatile charge pattern recognition globular C1Q domain in the C-terminus, engaging a wide range of ligands and trigger a serial of immune responses (Wang et al., 2013a;Song et al., 2015). Hemocyanins are multifunctional proteins existing in invertebrate hemolymph, which contribute to mollusks' innate immunity through phenoloxidase-like activities (Zhuang et al., 2015). Hence, DY may have transgressive performances in the complement system and hemocyanin-related functions. Interestingly, in aquatic species, unlike growth-related genes, which usually show overdominant expression patterns in the hybrids, the immunity-related genes seem more likely to be dominant or partial-dominant expressed in the hybrids (Yan et al., 2017;Zhang et al., 2017;Yang et al., 2018;Shen et al., 2020). This phenomenon might be due to aquatic hybrids' immune system characteristic, which needs further investigation.
(3) Stress response The stress response molecules, including antioxidant enzymes and acute-phase molecules, protect organisms' cells from the damage of excess reactive oxygen species (ROS) or superoxide anion (O 2− ) generating upon environmental stressors, or help organisms maintain regular molecular and cellular homeostasis (Aguirre et al., 2005;Vaák and Meloni, 2011;Wang et al., 2013b). In our study, the "environmental information processing" pathways were enriched from the ELD-J and ELOD genes ( Figure 5). Environmental information processing is usually the first step of organisms' stress responses. These results suggest DY is JP-dominant or overdominant in stress response related functions. This hypothesis is further supported by the expression of the 29 stress response genes in the hybrid. DY showed JP-dominant or overdominant expression in most gene groups. Notably, the CAT, SEGPX, and MGST genes were highly expressed in DY compared with the MPV. These genes are antioxidant enzyme genes whose functions include removing excess ROS (Chelikani et al., 2004), reducing organic hydroperoxide and hydrogen peroxidase (Arthur, 2000) , and catalyzing the nucleophilic attack of the glutathione's sulfur atom in organisms (Blanchette et al., 2007). Oxidative stress is common for mollusks as they mostly live in highvelocity water environments. Hence, the ability to reduce oxidative stress is crucial for mollusks to maintain physiological equilibrium. Similarly, in the abalone hybrid of H. discus hannai and H. gigantea, three antioxidant enzymes were found to be overdominant expressed (Di et al., 2015). The positive relationship between the hybrid abalone's stronger antioxidant functions and its superior survival rates could be concluded.

CONCLUSION
This study investigated the transcriptomes and miRNAomes of two H. diversicolor geographic genotypes and their hybrid (Figure 11). Evident non-additivity, ELD, and ELOD expression patterns were identified in the DEGs and DEMs, showing the widespread existence of dominant and overdominant models in the hybrid's heterosis formation. Functional analysis of the ELD and ELOD genes illustrated DY had maternal growth-related functions and paternal stress-response functions. Meanwhile, DY was overdominant in several functional pathways. The existence of miRNA regulation was found in several important functional pathways. The role of FGF, C1Q, HC, CAT, SEGPX, and MGST genes in the abalone heterosis generation is noteworthy. Our study provides a case of molecular heterosis study in mollusk species, which may help disclose the mysteries underlying this complex phenomenon.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA database under BioProject accession PRJNA701659.

AUTHOR CONTRIBUTIONS
SL conducted the experiments and wrote the manuscript. WY helped in the abalone culture, experiments design, and manuscript revision. XL helped in the abalone culture and experiments design. JK helped in the abalone culture. MH helped in the lab experiments. YG helped in the manuscript revision. CK sponsored the study and participated in the experiments design. All authors contributed to the article and approved the submitted version.