Transcriptome analysis of growth heterosis in pearl oyster Pinctada fucata martensii

Heterosis improves growth and survival of shellfish species. Although breeders have widely exploited heterosis, its underlying molecular mechanisms remain unclear. In this study, a 2 × 2 complete diallel cross was facilitated between two full‐sib families to produce two inbred families (A and D) and their reciprocal hybrid families (B and C) of pearl oyster Pinctada fucata martensii. Growth traits of the four families were compared at the adult stages. Transcriptome analysis was conducted on the four families using an Illumina sequencing platform. The results revealed that the growth traits of the four families significantly varied (P < 0.05). The mid‐parent heterosis values of shell length, shell height, shell width, shell weight, and total weight were 12.9%, 14.9%, 18.2%, 17.2%, and 33.2%, respectively. The B‐ and C‐inbred (A and D) triads had 79 and 68 differentially expressed genes (DEGs), respectively, which were dominantly nonadditive, including overdominance, underdominance, and low‐parent dominance. Gene ontology term analysis showed that the DEGs in the B‐ and C‐inbred triads were enriched for metabolic process, cellular process cell part, binding, and catalytic activity. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis indicated that the DEGs in the B‐ and C‐inbred triads were involved in focal adhesion, the P13K‐Akt signaling pathway, the mRNA surveillance pathway, and the focal adhesion pathway. The reliability of the sequencing data was confirmed by real‐time polymerase chain reaction analysis of six growth‐related genes. The findings of this study provide new insights into heterosis for growth traits and the design of genetic breeding programs for this species.

Heterosis improves growth and survival of shellfish species. Although breeders have widely exploited heterosis, its underlying molecular mechanisms remain unclear. In this study, a 2 9 2 complete diallel cross was facilitated between two full-sib families to produce two inbred families (A and D) and their reciprocal hybrid families (B and C) of pearl oyster Pinctada fucata martensii. Growth traits of the four families were compared at the adult stages. Transcriptome analysis was conducted on the four families using an Illumina sequencing platform. The results revealed that the growth traits of the four families significantly varied (P < 0.05). The midparent heterosis values of shell length, shell height, shell width, shell weight, and total weight were 12.9%, 14.9%, 18.2%, 17.2%, and 33.2%, respectively. The B-and C-inbred (A and D) triads had 79 and 68 differentially expressed genes (DEGs), respectively, which were dominantly nonadditive, including overdominance, underdominance, and low-parent dominance. Gene ontology term analysis showed that the DEGs in the B-and C-inbred triads were enriched for metabolic process, cellular process cell part, binding, and catalytic activity. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis indicated that the DEGs in the B-and C-inbred triads were involved in focal adhesion, the P13K-Akt signaling pathway, the mRNA surveillance pathway, and the focal adhesion pathway. The reliability of the sequencing data was confirmed by real-time polymerase chain reaction analysis of six growth-related genes. The findings of this study provide new insights into heterosis for growth traits and the design of genetic breeding programs for this species.
Pearl oyster Pinctada fucata martensii is one of the most important components of molluscan mariculture in southern China, particularly in Guangdong, Guangxi, and Hainan provinces. This species is cultured to produce nucleated round pearls referred to as the 'South China Sea Pearl'. However, pearl oyster cultures in these provinces have recently suffered from severe slow growth and mass mortalities. Breeding studies have been initiated to improve growth and resistance to diseases and increasingly deteriorating environments [1]. These studies include mass selection [1][2][3] and population hybridization [4]. Mass selection is an effective approach to improve the growth traits of P. fucata martensii. The fifth generation line selected for the faster growth rate of P. fucata martensii displayed a larger mean shell length and shell height than those of the control group [5]. Hybridization is an alternative approach to improve the growth traits of P. fucata martensii by crossing geographically isolated populations and selected lines [4].
Heterosis or hybrid vigor refers to the phenomenon in which a hybrid offspring exhibits phenotypic superiority over its parents. It has been commonly utilized to improve the growth traits and survival of shellfish species [4,[6][7][8]. Therefore, the genetic mechanisms of heterosis have gained considerable interest in genetics and breeding research. Although heterosis has been widely exploited by crop and animal breeders, its underlying molecular mechanisms remain poorly understood. Three classical hypotheses have been proposed to explain heterosis: dominance, overdominance, and epistasis. Dominance refers to the complementation by superior dominant alleles in a heterozygous hybrid. Overdominance ascribes allelic interactions at one or more loci for heterotic traits. Epistasis attributes heterosis to the epistatic interaction of beneficial nonallelic genes at two or more loci in hybrids [9].
With the rapid improvements and advancements in molecular biological methods, heterosis in crop and animals has been elucidated with molecular evidence. At present, several techniques are available for various aquatic species. These techniques include suppression subtractive hybridization, mRNA differential display techniques, and cDNA-amplified fragment length polymorphisms. A high-throughput sequencing of transcriptomes has been recently used to explore the molecular mechanism underlying heterosis [9]. Highthroughput sequencing of the transcriptome analyses of P. fucata martensii identified numerous candidate genes involved in embryogenesis, growth traits, and biomerialization [10,11]. Recent works have reported on the genomes of P. fucata martensii, which served as basis for investigating of their important traits [11].
We have recently developed a pearl oyster breeding program to improve the growth performance of cultured stocks, and a large number of families have been developed. In a previous study, we conducted an experiment to compare several gene expression profiles between purebred and hybrid families by qRT-PCR and found significantly different gene expression levels among the families [12]. In this study, the whole transcriptomes of inbred and hybrid families were sequenced to analyze the functional differences among the four families and provide a basis for investigating the mechanisms of growth heterosis in the species.

Experimental animals
A 2 9 2 complete diallel cross was made between two fullsib families (M and N) to produce four families. The four families were named as follows: (a) A, an inbred family produced by sister-brother mating from family M; (b) B, a hybrid family with the female parent from family M and the male parent from family N; (c) C, a reciprocal hybrid family with the female parent from family N and the male parent from family M; and (d) D, an inbred family produced by sister-brother mating from family N.
The inbred and hybrid families were developed in March 2014. In the complete diallel cross, a female was mated to a full-sib brother and an unrelated male to produce inbred and hybrid families. The larvae were reared following the techniques of Deng et al. [1]. The larvae were reared in 1000-L polyethylene tanks. The density was maintained at 1 individualÁmL À1 . Water temperature was 27.2 AE 1.2°C, and salinity was 30.0 AE 0.5 ppt. Daily feeding consisted of Isochrysis zhanjiangensis from days 2 to 7 and a mixture of I. zhanjiangensis and Chlorella sp. from days 7 to 55. On Day 21, plastic plates were provided as substrate for metamorphosis. On Day 50, individuals as large as 3-5 mm were removed from the plates, transferred to net-cages, and suspended in the sea of Liushagang, Zhanjiang. The shells were cleaned and placed in new nets at appropriate intervals.

Growth measurement and sampling
A total of 100 individuals with 2 years of age were sampled from each family. Shell length, shell height, shell width, and shell weight of each individual were measured. Ten individuals were sampled from each family, and adductor muscle of each individual was dissected. Each sample was formed by pooling adductor muscles of equal size from two individuals within a family and then stored at À80°C.

Total RNA extraction, cDNA library construction
Total RNA of the sample was extracted with TRIzol reagent. Trichloromethane was used for protein denaturation and phenol extraction. Isopropanol was used to precipitate the nucleic acids, which were then dissolved in DEPC water. RNase-free tips and tubes were used in RNA extraction at all times. The concentration of each sample as well as OD260/OD280 and OD260/OD230 was measured using SimpliNan. The total RNA, including RNA concentration, RIN value, 28S/18S, fragment length distribution, and cDNA library construction, was detected by an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit). Two samples were used for each family, and eight cDNA libraries were constructed.

Data analysis of RNA-seq
The project used internal software SOAPnuke to filter reads and remove reads with adaptors. Reads with unknown bases (N) were more than 5%. Low-quality reads, which were defined as a read in which the percentage of the base is less than 15, were greater than 20%. Clean data were obtained in FASTQ format. The statistics included Q20, Q30, and the ratio of clean reads. The clean reads were used in the subsequent analysis. The clean reads were mapped to the reference genome using HISAT after the reads were filtered. The bam format date of the reads that contrasted with the genome sequences was obtained [11]. The bam format date can be used by IGV. IGV can be compared with multiple samples to find the distribution of reads in exons, introns, UTRs, and intergenic regions.

Differential expression analysis
Differentially expressed genes (DEGs) were detected by the DEGSeq method. [log 2 Ratio (BP/SP)] > 1 and FDR≤0.001 were used as the threshold for a significant differential expression. On the basis of the DEGs, the classification and aggregation analysis of gene ontology (GO) function were performed and were divided into three major functions: cellular component, molecular function, and biological process. KEGG is a database resource from molecular level information to understand advanced functions and biological systems [13].
According to the method of Rapp et al. [14], the DEGs in a hybrid-inbred (parent) triad were classified into five possible expression classes of differential expression, that is, the expression levels of genes equal to the mid-parent (additivity), the high-or low-parent (high-or low-parent dominance), above the high-parent (overdominance), or below the low-parent (underdominance) [14]. The mid-parent value was calculated for two replicates and compared with the average hybrid expression for the replicates. A two-tailed homoscedastic t-test was performed. Genes with P > 0.05 were considered to be additively expressed, whereas genes with P < 0.05 were regarded as nonadditively expressed.

Quantitative real-time PCR
Gene expression profiles were validated by quantitative realtime PCR (qRT-PCR). The gene-specific primers for qRT-PCR were designed using PRIMER 5.0 (Premier Biosoft International, Palo Alto, CA, USA) ( Table 1). Ten individuals were used for each of the four families. The total RNA was extracted using TRIzol (Thermo Scientific TM). cDNA was synthesized using PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) [15,16]. A 5 lL PCR mixture consisted of 2.5 lL of SYBR mixture (TaKaRa), 0.2 lL of cDNA, 0.2 lL of each PCR primer, and 2.1 lL of ddH 2 O. The PCR amplification procedure was 40 cycles of 50°C for 2 min, 95°C for 2 min, 95°C for 15 s, 60°C for 1 min, 95°C for 15 s, 60°C for 1 min, 95°C for 30 s, and 60°C for 15 s. The results were expressed by 2 ÀMMCt . Each sample consisted of three replicates. The relative expression of the genes was calculated by the 2 ÀMMCt method, and the standard deviation was calculated among three replicates [17].

Statistical analysis
One-way ANOVA followed by Tukey test was performed to compare the differences in the mean values of growth traits and gene expression of the four families. Significance level was set at P < 0.05.
Mid-parent heterosis (MPH) was calculated for shell length, shell height, shell width, shell weight, and total weight of the samples using the following equation: MPH(%) = [(F 1 ÀMP) 9 100]/MP, where F 1 is the mean shell length, shell height, shell width, shell weight, or total weight of two reciprocal hybrids; MP = (P 1 + P 2 )/2, where P 1 and P 2 are the mean shell length, shell height, shell width, shell weight, and total weight of the two inbred families.

Comparison of the growth traits of the four families
Significant differences were detected in the measured growth traits of the four families (P < 0.05). These differences were observed between the inbred families, between the hybrid families, and between the inbred and hybrid families. Among the four families, hybrid C had the highest values for shell length, shell height, shell width, shell weight, and total weight. The MPH varied for the different growth traits with the values of shell length, shell height, shell width, shell weight, and total weight were 12.9%, 14.9%, 18.2%, 17.2%, and 33.2%, respectively (Table 2).

Transcriptome profiling of two hybrids and their parental families
After filtering low-quality, joint contamination, and high-content reads with an unknown base N, clean data were mapped to the genome using HISAT. The average ratios of alignment and unique mapping were 73.86% and 46.65% (Table 3). A total of 32 526 genes were detected, which included 28 156 known genes and 4430 predicted new genes. A total of 46 370 new transcripts were detected, of which 21 664 new alternative splice variants belonged to protein-coding genes, 4430 belonged to new protein-coding genes, and the remaining 20 276 belonged to long-chain noncoding RNAs.

DEGs between inbred and hybrid families
We investigated the variation in DEGs by performing pairwise comparisons between the inbred families and their reciprocal hybrid families. Differences were observed in gene expression between the two inbred families (A and D). A total of 594 DEGs (346 were upregulated and 248 were downregulated) were identified between A and D. Meanwhile, 283 DEGs (47 upregulated and 236 downregulated) were identified between reciprocal hybrid families B and C (Table 4).
We also investigated the DEGs in the B-and Cinbred triads. In the B-inbred triad, 79 DEGs were identified, of which 42 DEGs were upregulated and 37 DEGs were downregulated. In the C-inbred triad, 68 DEGs were identified, of which 49 DEGs were upregulated and 19 DEGs were downregulated (Table 4).

Functional analysis of DEGs
The DEGs between the two inbred families (A and D) presented genetic differences in parental breeders. GO analysis showed that the DEGs between the two   inbred families were enriched in the cellular and metabolic processes in biological processes; in the membrane and organelle of the cellular component; and in the binding and catalytic activities in molecular function. The enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of the DEGs between the two inbred families included microRNAs in cancer, ECM-receptor interaction, PI3K-Akt signaling pathway, and so on (Table 5). Differentially expressed genes between the hybrid and inbred families may be involved in heterosis. GO analysis revealed that the DEGs in the B-inbred triad were involved in the cellular process and metabolic processes in biological processes; in the cell and cell part of the cellular component; and in the binding and catalytic activities in molecular function. The enriched KEGG pathway of the DEGs in the B-inbred triad indicated that the DEGs were dominantly enriched in the metabolic pathways, including the Focal adhesion and the PI3K-Akt signaling pathway (Fig. 2). The genes in the enriched KEGG pathway were involved in immune response (nine genes) and growth and development (13 genes).
Gene ontology analysis showed that the DEGs in the C-inbred triad were involved in the cellular process and metabolic processes in the biological processes; the cell part, organelle, and membrane in the cellular component; and the binding and catalytic activities in molecular function. The enriched KEGG pathway of the DEGs in the C-inbred triad included the Focal adhesion, PI3K-Akt signaling pathway, and mRNA surveillance pathway (Fig. 2). The enriched genes in the KEGG pathway were involved in growth and development (nine genes) and immune response (22 genes).

Validation of differentially expressed genes using qRT-PCR
Quantitative real-time PCR was utilized to verify the differential expression of genes identified by transcriptome analysis. Six growth-related genes were selected  for qRT-PCR analysis. These genes included insulinlike growth factor-binding protein (IGFBP), cell growth-regulating nucleolar protein (CGRNP), E3 ubiquitin-protein ligase MYCBP2 (MYCBP2), zinc protein (ZP), tubulin beta chain (TBC), and ribosomal RNA methyltransferase (rRm). The RT-qPCR expression patterns of the selected DEGs were consistent with the results of RNA-seq analysis, indicating the reliability and accuracy of the RNA-seq method used in this study ( Fig. 3 and Table 6).

MPH values of the inbred and hybrid families
Hybridization between geographically isolated populations, strains, or families is typically used to improve the growth performance and survival of cultured stocks in shellfish species, such as catarina scallop Argopecten circularis [18], and bay scallop Argopecten irradians irradians [19], Pacific oyster Crassostrea gigas [6,20], Pacific abalone Haliotis discus hannai [7,21], sea scallops Placopecten magellanicus [22], and clam Meretrix meretrix [8]. For example, Deng et al. [21] developed a complete diallel cross among three geographically isolated populations (Dalian and Qingdao of China and Miyagi of Japan) of Pacific abalone H. discus hannai. They reported that the MPH values of shell length, shell width, and total weight are in the ranges of 0.5% to 12.2%, 0.7% to 16.7%, and À1.7% to 44.8%, respectively [21]. By crossing three strains of C. gigas that were successively mass selected for two generation from three culture stocks collected from China, Japan, and South Korea, Kong et al. [20] reported that the MPH values of shell height, shell length, and the whole body weight are in the ranges of À6.07% to 22.35%, À10.46% to 29.74%, and À27.76% to 89.23%, respectively. In the present study, the MPH values (heterosis) of shell length, shell height, shell width, shell weight, and total weight were 12.9%, 14.9%, 18.2%, 19.3%, and 33.2%, respectively. These findings are consistent with those of previous studies on other shellfish species. The MPH values reported here were considerably higher than the results of Wang et al. [23]. They constructed a 2 9 2 complete diallel cross between Indian and Chinese Sanya wild populations of pearl oyster Pinctada martensii. The MPH values of shell height, shell length, shell width, and shell weight were 4.08%, 2.78%, 7.27%, and 6.81%, respectively [23]. The present results indicated that the heterosis effects in pearl oyster P. fucata martensii could be explored by family mating.

DEGs of the purebred and hybrid families
The genetic basis of heterosis is complex and has yet to be fully elucidated. Several studies have explored the mechanisms underlying heterosis via diallel cross design in shellfish species [4,6,8,21]. For example, in a 3 9 3 complete diallel cross by sampling parents from three full-sib families of the clam Meretrix meretrix, Dai et al. [8] reported that nonadditive genetic effects are the main causes of genetic variation for shell length and whole body weight. By designing a complete diallel cross between three geographically isolated populations (Dalian and Qingdao of China and Miyagi of Japan) of Pacific abalone H. discus hannai, Deng et al. [21] reported substantial variance in the specific combining ability for growth traits of Pacific abalone H. discus hannai. Hedgecock and Davis [6] designed four incomplete diallel crosses and factorial mating using a set of six to nine inbred lines as both male and female parents to analyze the yield heterosis in juvenile (seed) and harvest-sized adult oysters. In these experiments, the variance in yield was partitioned into additive and nonadditive genetic components, and the specific combining ability was found to be significant in all four crosses. These studies imply that nonadditive gene actions determine these growth traits of the shellfish species mentioned above. With the development of modern biological technologies, recent works have attempted to explain heterosis at the molecular level. The genetic base of hybrid performance can be mainly or partly due to the complementarity of additive and nonadditive genetic effects [24]. Several studies have analyzed heterosis- associated gene expression in crop and animals by comparing the expression patterns of selected genes in inbred lines and hybrids or by performing highthroughput gene expression analyses via microarray profiling and transcriptome analyses [25][26][27][28]. These studies reported that nonadditive gene expression is prevalent between purebred and hybrid families. Hedgecock et al. [28] investigated the gene expression patterns underlying the growth heterosis in two partially inbred (inbreeding coefficient = 0.375) and two hybrid larval populations of the Pacific oyster C. gigas that were produced by a reciprocal cross between two inbred families. They categorized 70% of genotypedependent patterns of gene expression into classical modes of gene action and found nonadditive patterns accounted for 98% of gene expression [28]. In the  present study, we found that nonadditive gene expression patterns were prevalent between hybrid B or C and its parents. Combined functional analysis of nonadditive genes and expression pattern analysis of genes enriched in GO terms indicated that nonadditive genes participated in biological process regulation, stimulus response, and transcriptional regulation activity. Consequently, the heterosis in the growth traits observed in this study might be induced by the differential expression of the genes in the hybrids.

Growth-and immune-related genes
Several studies have confirmed that a positive relationship exists between growth traits and pearl traits (pearl thickness and pearl weight) [29,30]. Pearl traits were unexpectedly enhanced by improving the growth traits of pearl oyster stocks. Thus, growth traits are fundamental factors that must be considered in pearl oyster breeding programs. Growth traits are quantitative traits that are determined by multiple genes. The mining of growth-related genes can provide insights into the underlying mechanisms of growth traits and assist in breeding design. Transcriptome sequencing has been proven efficient for mining growth-related genes in pearl oyster [10,11]. In this study, we found several growthrelated genes in the DEGs in hybrid B-or C-inbred triads, such as IGFs, GHSR, and zinc finger protein.
The IGFs are a type of multifunctional cell proliferation regulators that play key roles in regulating invertebrate metabolism, growth, and reproduction [31,32]. As IGFs binding proteins, IGFBPs maintain IGF in the circulation, mediate the regulation of IGF independent activities, and regulate the growth and development of bone and other tissues. Zhang et al. [33] found that IGFBPs are present in pearl oyster P. martensii and are involved in IPL endocrine system regulation by a feeding experiment. GHSR is a G protein-coupled receptor, the natural ligand of Ghrelin hormones, and is mainly expressed in the hypothalamus and the pituitary gland [34]. The GHSR system is an important pathway for regulating growth hormones, which can promote growth and development as well as immune digestion. Shuto et al. [35] demonstrated that GHSR regulates growth hormone secretion and food acquisition in the arcuate nucleus. Zinc finger protein is a transcription factor that is divided into nine types, namely, C2H2, C8, C6, C3HC4, C2HC5, CCCH, C2HC, C4HC3, and C4 [36][37][38]. CCCH zinc finger protein is involved in growth regulation and immune function. OMA-1 and OMA-2 in nematodes are involved in oocyte maturation and embryonic development [39,40]. In C. elegans, xc3 h-3b plays an important regulatory role in renal differentiation [41]. MCP-induced proteins, which comprise a new family of CCCH zinc finger proteins, regulate macrophage activation, reduce inflammation, and prevent diseases. In this study, we found that several growth-related genes were present in the DEGs in the hybrid B-and C-inbred triads, indicating that their gene expression contributed to the growth heterosis in pearl oyster P. fucata martensii.
Innate and adaptive immune responses are well-recognized host defense mechanisms of plants, fungi, and animals. Invertebrates rely solely on innate immune responses for defense against pathogens in natural environments. Identifying immune-related genes and their pathway can elucidate the mechanisms against pathogens and develop pathogen-resistant strains. In this study, we found several immune-related genes by analyzing DEGs in the hybrid B-and C-inbred triads, such as ubiquitin (Ub) protein ligase E3, histone H2A, and phytanoyl-CoA. Cheng et al. [42] identified a novel Ub protein ligase E3 (CgE3Rv1) from Crassostrea gigas. They found that CgE3Rv1 is a Ub protein ligase E3 that is involved in the immune response to LPS and in the interaction with cell surface signaling molecules of the neuroendocrine immune system in oysters. Pellino is a highly conserved E3 class Ub ligase that plays a role in the pathogenesis of WSSV and the antiviral mechanism of shrimp Litopenaeus vannamei [43]. Histones H2A, H2B, H3, and H4 were found in the shrimp Litopenaeus vannamei and catfish Parasilurus asotus, hematopoietic histone proteins displaying antimicrobial activity and Cathepsin D and metalloproteinases are involved in immune response [44,45]. A yeast two-hybrid system was used to search for potential proteins of FKBP52. The results indicated that PAHX is an important candidate for studying cell signaling pathways involving FKBP52 in the presence of immunosuppressant drugs [46]. The PhyH-like expression of Xenopus tadpoles was expressed in autoreactive immune cells that were distributed throughout the body during the refractory period [47].

Conclusion
The transcriptome analysis of hybrid families and their parental families were conducted using the RNA-seq technology. We found that DEGs in the hybrids and their parents might be associated with growth heterosis. Analyzing the gene action models between the two hybrids and their parents, we proposed overdominance might play important roles on growth heterosis. The transcriptome data from this study will help to understand of molecular mechanism underlying growth heterosis in pearl oyster P. fucada martensii.