- Email: [email protected]
Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂
Accepted Manuscript Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂
Guodong Zheng, Chenbin Wu, Juan Liu, Jie Chen, Shuming Zou PII: DOI: Article Number: Reference:
S0044-8486(19)30779-3 https://doi.org/10.1016/j.aquaculture.2019.734317 734317 AQUA 734317
To appear in:
aquaculture
Received date: Revised date: Accepted date:
2 April 2019 30 June 2019 16 July 2019
Please cite this article as: G. Zheng, C. Wu, J. Liu, et al., Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂, aquaculture, https://doi.org/10.1016/j.aquaculture.2019.734317
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂
PT
Guodong Zheng, Chenbin Wu, Juan Liu, Jie Chen, Shuming Zou*
RI
Key Laboratory of Freshwater Aquatic Genetic Resources, Ministry of Agriculture,
SC
Genetics and Breeding Center for Blunt Snout Bream, National Demonstration Center
NU
for Experimental Fisheries Science Education, Shanghai Ocean University, Shanghai 201306, China
MA
*Correspondence Authors: Dr. Shu-Ming Zou
Tel.: +86-021-61900345
ED
Fax: +86-021-61900345
AC C
EP T
Email: [email protected]
1
ACCEPTED MANUSCRIPT
Abstract A novel backcross variety (BC) [Megalobrama amblycephala (MA) ♀ × (Megalobrama amblycephala ♀ × Culter alburnus (CA) ♂) ♂] was generated by artificial insemination. The BC exhibits growth superiority and enhanced digestive
PT
enzyme activity than that of its parents. We conducted the transcriptome analysis
RI
using hepatopancreas between the BC and its parents. Three hundred and twenty-five
SC
and 872 differentially expressed genes (DEGs) were identified in ‘BC vs. MA’ and
NU
‘BC vs. CA’, respectively. The GO and KEGG enrichment analysis of DEGs were classified. From the 134 overlapping DEGs and their pathway, the digestive enzyme
MA
genes, such as TRY, ELA1, CTRL1, CPA2, and BAL, and the IGF system genes, such as IGF1, IGF2a, and IGFBP2b, of the BC had significant difference in expression.
ED
Furthermore, the expression of downstream genes related to the protein and fatty acid
EP T
synthesis pathways, such as PI3KR, RAPTOR, EIF4E, CS, MDH, FAS, ELOVL1, ELOVL5, and ELOVL6, were significantly up-regulated. In addition, a total of 40732
AC C
SSRs and 52951 SNPs were identified from the transcriptome. The results revealed that the enhanced digestive enzymes activity and synthesis of proteins/fatty acids might be contributed to growth superiority of backcross BC. Key words: transcriptome; backcross; growth superiority; digestive enzyme; protein/fatty acid synthesis
2
ACCEPTED MANUSCRIPT
1. Introduction Blunt snout bream (Megalobrama amblycephala, MA) and topmouth culter (Culter alburnus, CA) are freshwater fish species of the subfamily Cultrinae. In China, both species are important in the fish polyculture system. Although the M.
PT
amblycephala grows faster than the C. alburnus, but its meat quality is inferior to that
RI
of the C. alburnus. To produce new varieties with growth superiority and high quality
SC
meat, hybridization has been introduced into the artificial breeding of these two
NU
species. According to statistics, the yield of the hybrids between M. amblycephala and C. alburnus has exceeded 50 thousand tons in China in 2018 (China Aquaculture
MA
Network, www.shuichan.cc). In previous study, a fertlie hybrid lines between M. amblycephala and C. alburnus were established (Xiao et al., 2014), this hybrid were
ED
also found to have obvious heterosis in growth performance, feed utilization and
EP T
muscle quality (Zheng et al., 2015; Zheng et al., 2019). Moreover, the total amino acid and essential amino acid content of hybrids between M. amblycephala and C.
AC C
alburnus were higher than those of their parents (He et al., 2014), SRAP marker analysis showed that the genetic diversity of this hybrids increased significantly (Jia et al., 2011) and the number of intermuscular bones decreased significantly (Jiang et al., 2016). In addition, a similar hybrid between M. terminalis and C. alburnus also showed growth advantage (Guo et al., 2018). RNA-Seq has been widely used as an effective tool for transcriptome analysis in order to discover, profile, and quantify the RNA transcripts (Wang et al., 2009).
3
ACCEPTED MANUSCRIPT
RNA-Seq combines the advantages of high throughput, low background noise and high sensitivity, making it feasible to map transcribed regions, quantify gene expression, and distinguish different isoforms and allelic expression (Wang et al., 2009; Ozsolak and Milos, 2011). Recently, the differently expressed genes (DEGs)
PT
and pathways related to heterosis have been identified through the transcriptome
RI
profiles analysis in some commercially crops and fishes, such as hybrid maize (Zea
SC
mays L.) (Bi et al., 2014), hybrid rice (Oryza sativa) (Zhai et al., 2013), hybrid
NU
soybean (Glycine max (L.) Merr) (Zhang et al., 2017), hybrid pufferfish (Takifugu) (Gao et al., 2013), hybrid grouper (Epinephelus spp.) (Sun et al., 2016a, b), etc. These
MA
studies could help us to further understand the mechanisms underlying the heterosis in crops and fishes.
ED
In order to reveal the molecular mechanism of the heterosis of M. amblycephala
EP T
and C. alburnus, transcriptome analysis was carried out for the two species and their hybrids. It is hoped that effective DEGs can be screened to reasonably explain the
AC C
heterosis of fish, and a large number of molecular markers can be developed for the study of QTL analysis and genetic maps.
2. Materials and methods 2.1. Experimental fish and feeding trial M. amblycephala (MA) and C. alburnus (CA) broodstocks were obtained from the Genetics and Breeding Center for Blunt Snout Bream of Ministry of Agriculture,
4
ACCEPTED MANUSCRIPT Shanghai, China. MC hybrids (MA♀ × CA♂) were produced first, when MC hybrids reached sexual maturity, the backcross BC (MA♀ × MC♂) was produced and parental selfbred strains of MA (♀ × ♂) and CA (♀ × ♂) were also produced by artificial insemination. All fish used in the study came from a single cross family, eggs from
PT
the same female used for the backcross with the male hybrid also used for maintaining
RI
the parental strains. These three crosses (BC, MA and CA) with the same age and
SC
similar sizes were randomly distributed into 3 concrete ponds (6 × 4 × 1.2 m, L : W :
NU
H) with water depth of 1m. Each pond was stocked all three crosses and each cross had 50 fry (total of 150 fish per pond). The feeding period was 90 days. During the
MA
feeding trial, fish were fed twice (06:30 and 17:30) each day with equal amounts of floating compound feed (containing 33% protein, 3% lipid and 8.5% fiber, Tech-bank
ED
Co., Ltd, Ningbo, China) and until the feed is eaten completely by the fish. One half
EP T
of the pond water was changed weekly to ensure fresh water quality. Temperature of
AC C
water ranged from 22~30 ℃ and dissolved oxygen was above 5.0 mg/L.
2.2. Growth measurement and sample collection At the end of feeding trial, after 24 hours of starvation, the weight gain rate (WGR), specific growth rate (SGR), hepatosomatic index (HSI) and condition factor (CF) were calculated. Ten individuals per cross were euthanized with 100 mg/L of MS222 (Sigma, USA) before tissues were collected. Hepatopancreases (for RNA extraction, sequencing and digestive enzyme activity assays) was collected and stored
5
ACCEPTED MANUSCRIPT
immediately at -80°C.
2.3. Digestive enzyme activity assays The whole hepatopancreas of three samples per cross was weighed and
PT
homogenized (dilution 1:10) in ice-cold buffer. The homogenization buffer solution
RI
was 0.02 M Tris/0.01 M phosphate, pH 7.0, in 50% (v/v) glycerol (Moro et al., 2010).
SC
The extract was centrifuged at 3000 rpm for 10 min at 4°C, and the supernatant was
NU
used as the enzyme source using BSA (Sigma, USA) as standard. The activity of total protease was detected using Folin-phenol method and the activity of total lipase was
MA
detected using commercial kits (ref. no. A054) produced by Jiancheng Bioengineering
ED
Institute (Nanjing, China) (Li et al., 2012).
EP T
2.4. RNA extraction, library construction and sequencing Total RNA from hepatopancreas samples was extracted using TRIzol reagent
AC C
(Invitrogen, Carlsbad, CA, USA) and purified using the RNeasy Mini Kit (Qiagen, Valencia, CA, USA). RNA degradation and contamination were detected by 1% agarose gels. The integrity and concentration of RNA samples were examined by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Equal amounts of nine RNA samples (three samples per cross) were used for cDNA synthesis and RNA-seq. Before the construction of library, Magnetic Oligo-dT beads (Invitrogen, USA) were used to isolate poly (A) + mRNA, and then cDNA
6
ACCEPTED MANUSCRIPT libraries were constructed and sequenced by Illumina Hiseq™ 4000 platform (Majorbio Biotech Co., Ltd, Shanghai, China). Briefly, 8 μg of total RNA for each sample was used to construct libraries by using Truseq TM RNA sample prep Kit (Illumina, USA) according to the manufacturer’s protocol. The constructed DNA
PT
template was enriched by PCR amplification (15 cycles) . Amplicons were collected
RI
and purified by Certified Low Range Ultra Agarose (Bio-Rad, USA) gel
SC
electrophoresis. Before sequencing, the DNA libraries were quantified by using
NU
TBS-380 micro fluorometer with Picogreen® reagent (Invitrogen, USA). Clone clusters were generated on Illumina cBot, using Truseq PE Cluster Kit v3-cBot-HS,
MA
and high- throughput sequencing was performed on an Illumina Hiseq™ 4000
ED
sequencer, using Truseq SBS Kit (300 cycles).
EP T
2.5. De novo assembly and functional annotation of unigenes The raw paired end reads from all transcriptomes were cleaned by removing
AC C
adapter contamination, low quality sequences (reads with over 10% unknown base pairs‘N’) and empty reads. All clean reads were then assembled using the denovo assembly program Trinity. The assembled unigenes of three crosses were used for BlastX search and annotation against the NR, Swissprot, String (Franceschini et al., 2013), COG (Clusters of Orthologous Groups of proteins), KEGG (Kyoto Encyclopedia of Genes and Genomes) databases with an E-value cut-off of 10−5 . Functional annotation including biological process, molecular function and cellular
7
ACCEPTED MANUSCRIPT
component for gene ontology (GO) terms was analysed by using BLAST2GO software (Conesa et al., 2005).
2.6. Comparative expression analysis
PT
The read counts were further normalized into FPKM (Fragments Per Kilobase of
RI
exon model per Million mapped reads) values. The FPKM values from the three
SC
libraries were pairwise compared, and the fold changes (FC) of the genes were
NU
calculated by using RSEM software (version 1.2.7) (Li and Dewey, 2011). If |log2(FC)| ≥ 1 and false discovery rate (FDR) ≤ 0.05, these genes were considered as
MA
having significant differential expression (Anders and Huber, 2010). The expressions of unigenes were displayed intuitively by Volcano Plot. After DEGs in the
ED
hepatopancreas between backcrosses and their parents were screened, the enrichment
EP T
analysis of GO and KEGG pathways was carried out by Goatools software (version0.4.7) and KOBAS 2.0 software (Xie et al., 2011). The R program was used
AC C
to depict the heat maps for gene clustering in the present study.
2.7. Detection of SSRs and SNPs Assembled
transcriptomes
were
screened
for
microsatellites
using
Msatcommander program (Faircloth, 2008). The parameters were designed to identify perfect mono-hexa- nucleotide motifs with a minimum of ten repetitions for mononucleotides, six repetitions for dinucleotides, and five repetitions for tri- hexa
8
ACCEPTED MANUSCRIPT nucleotides. Potential SNPs were filtered using the program Samtools (Li, 2011) and VarScan v.2.2.7 (Koboldt et al., 2009 ) with the default parameters.
2.8. Qualitative real-time PCR (qRT-PCR)
PT
Total RNA of hepatopancreas was extracted with the method described above, and
RI
cDNA was synthesized by PrimeScript RT reagent Kit (Takara, Janpan). The cDNA
SC
was then used for qRT-PCR analysis using specific primers (Table 1). qRT-PCR was
NU
performed by using SYBR Green RT-PCR kit (Takara, Janpan) on the CFX96 Touch™ real- time PCR Detection System (BioRad, USA). The housekeeping gene
MA
β-actin had been proved to be stable between the comparison groups and then used as the internal control. All reactions were performed in triplicates. The relative
ED
expression ratio of the target genes versus β-actin gene was calculated using 2-ΔΔCT
EP T
method (Livak and Schmittgen, 2001).
AC C
2.9 Statistical analysis
Digestive enzyme activity and gene expression data were carried out by using SPSS version 17 (Michigan Avenue, Chicago, IL, USA) for significant differences. If significant (P < 0.05) differences were found in factors, Duncan's multiple range test (Duncan, 1955) was used to rank the means.
3. Result
9
ACCEPTED MANUSCRIPT
3.1. Growth performance and digestive enzyme activity The growth performances of the backcross BC, MA, and CA are presented in Table 2. The weight gain rate (WGR) of BC, MA, and CA was 680.00%, 589.38%, and 334.61%, respectively. The WGR of BC was 15.38% and 103.22% higher than
PT
that of its parents MA and CA, this result showed that the growth superiority of the
RI
BC was significant (P < 0.05). The condition factor (CF) of BC was slightly and
SC
significantly (P < 0.05) higher than that of MA and CA, respectively. However,
NU
hepatosomatic index (HSI) of BC was intermediated between their two parents. In addition, the protease and lipase activity in hepatopancreas of BC was significantly (P
MA
< 0.05) higher than that of their parents, protease activity of CA was slightly higher
ED
than that of MA and lipase activity of MA was higher than that of CA (P < 0.05).
EP T
3.2. High-throughput sequencing and mapping to the reference genome Over 7.45 Gb of raw data for each sample and 159 million 100-bp paired-end
AC C
reads with an average of 53 million reads for each of the three samples (Table 3) were obtained. The raw data were cleaned and quality checked, and then assembled. Approximately, 49 million (BC), 55 million (MA), and 51 million (CA) clean reads were produced; 83.28% (BC), 89.64% (MA), and 89.38% (CA) of these clean reads were mapped to the assembled sequences (≤ 5 base mismatches) (Table 3). An assembly of the reads generated 77102 (BC), 63367 (MA) and 67087 (CA) transcripts with an average length of 923 bp, 835 bp and 842 bp, respectively. And a
10
ACCEPTED MANUSCRIPT
total of 62248 (BC), 55182 (MA) and 57366 (CA) unigenes were generated, with an average length of 813 bp, 745 bp and 747 bp, respectively (Table 4). An average of 61.78% of the transcripts and 66.29% of the unigenes were < 600 bp in size for the three samples (Supplementary Fig. 1). The results showed that ORFs of 41713
PT
unigenes (67.01%) in BC, 34654 unigenes (62.79%) in MA and 34602 unigenes
RI
(60.32%) in CA were successfully predicted using Trinity respectively. Overall, there
SC
were 25 classifications produced from the COG/KOG/NOG-annotated putative
NU
proteins, including ‘metabolism’, ‘cellular processes a nd signaling’ and ‘information storage and processing’, in accordance with the GO annotation categories
MA
(Supplementary Fig. 2 and Supplementary Table 1).
Characterization of BC unigenes by searching against public databases: E-value
ED
distribution of the top hits in the databases showed 11205 (35.39%) unigenes with a
EP T
strong homology (<1.0e-50) (Fig. 1A), 27629 (87.28%) unigenes had a similarity higher than 80%, while 36754 (11.61%) showed a similarity between 60% and 80%
AC C
with respect to the identity distribution pattern. Therefore, 98.89% of the unigenes showing an identity higher than 60% along with a high-quality similar distribution supported the reliability of the de novo assembly (Fig. 1B). A total of 31657 putative known unigenes were found in all reference species, 16903 (53.39%) were found in Danio rerio (Fig. 1C).
3.3. Identification and analysis of the diff erently expressed genes (DEGs)
11
ACCEPTED MANUSCRIPT
FPKMs of each gene in the hepatopancreas of the backcross BC were compared with those of their parents. The DEGs of ‘BC vs. MA’ and ‘BC vs. CA’ were identified by filtering the data with the following criteria: absolute log2 (FC) ≥ 1, FDR < 0.05, and FPKM > 2. According to these criteria, volcano plots were used to
PT
describe the number of significant up and downregulated DEGs of ‘BC vs. MA’ (Fig.
RI
2A) and ‘BC vs. CA’ (Fig. 2B).
SC
As shown in Fig. 3, Venn diagram showed that there were 325 DEGs (130
NU
up-regulated and 195 down-regulated) in ‘BC vs. MA’ and 872 DEGs in (388 up-regulated and 484 down-regulated) ‘BC vs. CA’. Among those genes, 134
MA
overlapping DEGs were found in ‘BC vs. MA’ and ‘BC vs. CA’. Based on GO annotation, the hepatopancreas DEGs of ‘BC vs. MA’ and ‘BC vs. CA’ were enriched
ED
into 46 and 53 GO terms respectively, which provides an overview of ontology
EP T
content. Metabolic process, cell part and binding were the most three enriched GO terms in the biological process, cellular component and molecular function categories,
AC C
respectively (Fig. 4A and B). By KEGG analysis, the DEGs of ‘BC vs. MA’ were classified to 5 categories with 54 pathway and DEGs of ‘BC vs. CA’ were classified to 5 categories with 51 pathway (Supplementary Fig. 3-4), in which the protein digestion-absorption pathway and pancreatic secretion pathway were the most significant pathway in both ‘BC vs. MA’ and ‘BC vs. CA’.
3.4. Analysis of diff erently expressed genes related to growth
12
ACCEPTED MANUSCRIPT
By analyzing these 134 overlapping DEGs and their pathway, we found that 5 genes, including trypsinogen (TRY), pancreatic elastase 1 (ELA1), chymotrypsin- like 1 (CTRL1), carboxypeptidase 2 (CPA2) and bile salt-activated lipase (BAL) of BC had higher FPKMs value than both MA and CA, which provides further evidence that
PT
digestive enzyme genes were more active in the BC (Fig. 5 and Table 5).
RI
From these 134 overlapping DEGs and their pathway, we also analyzed the
SC
insulin- like growth factor (IGF) system related genes in hepatopancreases, such as
NU
IGF1, IGF2a, IGF binding protein1 (IGFBP1) and IGFBP2b were differentially expressed in BC when compared with those of its parents (Fig. 5 and Table 5). In
including
MA
addition, the DEGs were identified in protein and fatty acid synthesis pathways, phosphoinositide-3-kinase
regulatory
subunit
(PI3KR),
regulatory
ED
associated protein of mTOR (RAPTOR), eukaryotic translation initiation factor 4E
EP T
(EIF4E), citrate synthase (CS), malate dehydrogenase (MDH), fatty acid synthetase (FAS) and elongation of very long chain fatty acids protein (ELOVL1, ELOVL5 and
AC C
ELOVL6) (Fig. 5 and Table 5).
3.5. SSRs and SNP discovery Using the Msatcommander program, a total of 12718, 12922 and 15092 potential SSRs contained repeats of more than one nucleotides with a minimum of five repetitions were identified from MA, CA and BC, respectively (Fig. 6). The number of SSRs with ten repetitions is the largest in all three crosses and the BC group has
13
ACCEPTED MANUSCRIPT
more SSRs than MA and CA group. Among these SSRs, GT/GT, GAT and GATA were most common in dinucleotide, trinucleotide and quadnucleotide SSRs, respectively. A total of 52951 candidate SNPs were identified in 3 crosses. Of these SNP
PT
candidates, 38787 SNPs were putative transitions (Ts), and 14164 SNPs were putative
RI
transversions (Tv) with a mean Ts: Tv ratio of 2.73 (Fig. 7). The SNPs were then
SC
categorized into 4 classes, including class 1 (C/A, A/C, T/G and G/T) for 11.44%,
NU
class 2 (C/T, G/A, T/C and A/G) for 73.25%, class3 (C/G and G/C) for 6.84%, and
MA
class 4 (A/T and T/A) for 8.47%.
3.6. Validation of the RNA-seq analysis by qRT-PCR
ED
To validate the DEGs in FPKM values, ten DEGs were randomly selected for
EP T
qRT-PCR assay in hepatopancreas of these three crosses. The results were further compared with the FPKM values generated from RNA-seq. Our results showed that
AC C
the data of qRT-PCR were consistent with those of RNA-seq (Fig. 8).
4. Discussion
Hybridization is widely used to improve growth rate or other economic traits of animals and plants (Lafarga-De la Cruz et al., 2013; Zheng et al., 2017; Chen et al., 2018; Wang et al., 2019; Zheng et al., 2019). In recent years, many studies have focused on exploring the molecular basis of heterosis in hybrid species by
14
ACCEPTED MANUSCRIPT
transcriptome analysis (Zhang et al., 2008; Gao et al., 2013; Zhai et al., 2013; Sun et al., 2016a, b). In the present study, comparative transcriptomic analysis was performed on hepatopancreas of the backcross BC and their parents MA and CA. Among the annotated transcripts, 1197 DEGs (325 and 872 DEGs in ‘BC vs. MA’ and
PT
‘BC vs. CA’, respectively) were identified, and 134 overlapping DEGs were
RI
confirmed. By analyzing these 134 overlapping DEGs and their pathway, many
SC
growth superiority-related genes were found to be diff erentially expressed in BC
NU
compared with MA and CA. Based on the function of these DEGs and their location in the pathway, a overview map of the generation of growth superiority was produced
MA
(Fig. 9).
ED
4.1. Gene expression and activity of digestive enzyme
EP T
The growth and development of fish depends on nutrients, especially amino acids and fatty acids. High-quality digestive juices can break down food efficiently,
AC C
improving the feed utilization of fish (Wang and Hartsuck, 1993; Jayaraman et al., 2006). In the present study, we found that the expression of four protein digestive enzyme genes, such as TRY, ELA1, CTRL1 and CPA2, and one fat digestive enzyme gene BAL was higher in the BC than that of their parents (Fig. 9). Furthermore, we found that the activity of most protease and lipase also increased significantly in BC, which was consistent with digestive enzyme genes expression above. Some studies showed that the growth performance of fish was related to digestive enzyme activity
15
ACCEPTED MANUSCRIPT
(Salze et al., 2012; Lu et al., 2015). With increased activity of digestive enzymes, large molecules of protein and fat can be broken down more efficiently into small molecules of amino acids and fatty acids, which can provide more raw materials for
PT
protein/fatty acid synthesis.
RI
4.2. IGF system and protein synthesis
SC
The growth performance of vertebrates or fish is primarily effected by the
NU
growing axis (Duan 1997; Duan 1998). In this axis, IGF system is involved in cell regeneration, proliferation and differentiation (Stewart and Rotwein 1996; Mahardini
MA
et al., 2018). In the present study, four DEGs associated with IGF system were found between BC and their parents (Fig. 9). Among them, the expression of IGF1 and
ED
IGF2a in the hepatopancreas of BC was significantly up-regulated than that of their
EP T
parents, their main function was to amplify and activate downstream pathways (such as protein synthesis pathway mentioned in next paragraph) (Sun et al., 2016a).
AC C
However, IGFBP1 exhibited a significantly decreased expression in the BC compared with their parents. This might indicate that reduced IGFBP1 protein allow IGFs to be released from the IGFBP-IGFs conjugates, therefore, more separated IGFs could combine with IGFR on cell surface and exert their role, which is similar to previous study in Scophthalmus maximu (Hu et al., 2012). On the contrary, the expression of the IGFBP2b mRNA in the BC was up-regulated significantly, the same results were also found in hybrid grouper (Sun et al., 2016a).
16
ACCEPTED MANUSCRIPT
In present study, besides the increased expression of IGF1 and IGF2a, the genes PI3KR, RAPTOR and EIF4E in the protein synthesis pathway also had higher expression in the BC, indicating the strengthened ability of protein synthesis in the BC. This results are consistent with previous study of hybrid grouper (Sun et al.,
PT
2016a). Combined with the above paragraph, many phosphorylation events, including
RI
PI3K/AKT pathway, are activated by the binding of IGF-1 to its receptor IGFR1 (Stitt
SC
et al., 2004), thus leading to a series of anabolic effects, such as protein synthesis
NU
(Burgos and Cant, 2010) and glycogen synthesis (KallooHosein et al., 1997). When stimulated by IGF-1, the PI3K/AKT signal activates mTOR (Bibollet-Bahena and
MA
Almazan, 2009), mTORC1 is then formed by the combination of Raptor with mTOR, which can phosphorylate EIF4E-binding protein 1 (EIF4E-BP1) (Ma and Blenis,
ED
2009). The phosphorylated EIF4E-BP1 can promote the dissociation of EIF4E from
EP T
the inhibitory complex of EIF4E and EIF4E-BP1 (Kimball et al., 1999). Finally, the
AC C
increased available EIF4E improves the synthesis of protein (Kimball et al., 1999).
4.3. Fatty acid synthesis The expression of many genes involved in the fatty acid synthesis pathway, such as CS, MDH, FAS, ELOVL1, ELOVL5 and ELOVL6 was up-regulated in the hepatopancreas of backcross BC compared with those in their parents (Fig. 9). These up-regulated genes may lead to faster rate of fatty acid synthesis and their specific roles in fatty acid synthesis pathway were explained as follows. As the direct raw
17
ACCEPTED MANUSCRIPT
material for fatty acid synthesis, acetyl CoA can bind with oxaloacetate to generate citric acid by the catalysis of citrate synthetase (CS), and enter the cytoplasm from mitochondrion (Verma et al., 2018; Xu et al., 2017). Citric acid will be cleaved into acetyl CoA and oxaloacetate in cytoplasm. Oxaloacetate then has to enter the
PT
mitochondrion by the catalysis of MDH for the next cycle of combining acetyl CoA.
RI
In cytoplasm, acetyl CoA can turn into malonyl Co A by catalysis of acetyl CoA
SC
carboxylase (ACC), fatty acid synthetase (FAS) then catalyzes the transformation of
NU
both acetyl CoA and malonyl CoA into fatty acid (16C) (Smith et al., 2003; Leng et al., 2012). Additionally, the ELOVLs conduct the extension of very long chain fatty
MA
acid (18C or 20C or 24C) in endoplasmic reticulum (ER) to play a variety of
EP T
4.4. SSRs and SNPs
ED
biological roles (Morais et al., 2009; Gregory et al., 2010; Green and Olson, 2011).
By transcriptome analysis, a large number of SSRs and SNPs have been identified
AC C
in BC and their parents. We found that SSRs with 10 repeats and SNPs of transition (Ts) were the most numerous, respectively (Li et al., 2018). However, the number of SSRs and SNPs in BC was significantly higher than that of their parents, this might be because the backcross BC genome fused the genomes of both parents. This SSRs and SNPs will provide a helpful resource for marker-assisted selection (MAS), quantitative trait locus (QTL) association and population genetics analysis (Li et al., 2018).
18
ACCEPTED MANUSCRIPT
5. Conclusion In the present study, we performed a hepatopancreas transcriptome analysis for a backcross BC and their parents to explore the molecular mechanisms of growth
PT
superiority. The results revealed that the digestive enzymes activity and the synthesis
RI
of proteins/fatty acids might be important factors for growth superiority of BC.
SC
Moreover, a large number of SSRs were identified in hepatopancreas transcriptome,
NU
which would be beneficial to QTL analysis and genetic linkage. Overall, these
Acknowledgments
ED
growth superiority of hybrid fish.
MA
findings provide new insights into molecular and physiological basis underlying the
EP T
This work was supported by grants from the National Natural Sc ience Foundation of China (31572220), Project funded by China Postdoctoral Science Foundation
AC C
(2019M651473), and the Shanghai University Knowledge Service Platform (ZF1206).
19
SC
RI
PT
ACCEPTED MANUSCRIPT
NU
References Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome
MA
Biol. 11, R106.
Bi, Y.M., Meyer, A., Downs, G.S., Shi, X.J., El-Kereamy, A., Lukens, L., Rothstein, S.J., 2014. High throughput RNA sequencing of a hybrid maize and its parents shows different mechanisms responsive to nitrogen limitation. Bmc Genomics 15, 77.
ED
Bibollet-Bahena, O., Almazan, G., 2009. IGF-1-stimulated protein synthesis in oligodendrocyte progenitors requires PI3K/mTOR/Akt and MEK/ERK pathways. J. Neurochem. 109,
EP T
1440-1451. Burgos, S.A., Cant, J.P., 2010. IGF-1 stimulates protein synthesis by enhanced signaling through mTORC1 in bovine mammary epithelial cells. Domest. Anim. Endocrin. 38, 211-221. Chen, J., Luo, M., Li, S.N., Tao, M., Ye, X.L., Duan, W., Zhang, C., Qin, Q.B., Xiao, J., Liu, S.J., 2018. A comparative study of distant hybridization in plants and animals. Sci. China Life Sci.
AC C
61, 285-309.
Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analys is in functional genomics research. Bioinformatics 21, 3674-3676. Duan, C.M., 1997. The insulin-like growth factor system and its biological actions in fish. Am. Zool. 37, 491-503. Duan, C.M., 1998. Nutritional and developmental regulation of insulin-like growth factors in fish. J. Nutr .128, 306s-314s. Faircloth, B.C., 2008. Msatcommander: detection of microsatellite repeat arrays and automated, locus-specific primer design. Mol. Ecol. Resour. 8, 92-94. Firth, S.M., Baxter, R.C., 2002. Cellular actions of the insulin-like growth factor binding proteins. Endocr. Rev. 23, 824-854. Franceschini, A., Szklarczyk, D., Frankild, S., Kuhn, M., Simonovic, M., Roth, A., Lin, J.Y.,
20
ACCEPTED MANUSCRIPT Minguez, P., Bork, P., Von Mering C., Jensen, L.J., 2013. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 41, D808-D815. Gao, Y., Zhang, H., Gao, Q., Wang, L., Zhang, F., Siva, V.S., Zhou, Z., Song, L., Zhang, S., 2013. Transcriptome analysis of artificial hybrid pufferfish Jiyan-1 and its parental species: implications for pufferfish heterosis. PLoS One 8, e58453. Green, C.D., Olson, L.K., 2011. Modulation of palmitate-induced endoplasmic reticulum stress and apoptosis in pancreatic beta-cells by stearoyl-CoA desaturase and Elovl6. Am. J. Physiol.
PT
Endocrinol. Metab. 300, E640-649. Gregory, M.K., See, V.H., Gibson, R.A., Schuller, K.A., 2010. Cloning and functional characterisation of a fatty acyl elongase from southern bluefin tuna (Thunnus maccoyii).
RI
Comp. Biochem. Phys. B 155, 178-185.
Guo, H.H., Zheng, G.D., Wu, C.B., Chen, J., Jiang, X.Y., Zou, S.M., 2018. Comparative analys is
SC
of the growth performance and intermuscular bone traits in F 1 hybrids of black bream (Megalobrama terminalis) (♀) × topmouth culter (Culter alburnus) (♂). Aquaculture 492, 15-23.
NU
He, Z.L., Liu, S.J., Xiao J., Hu, F.Z., Wen, M., Ye, L.H., Zhang, C., Xu, K., Tao, M., 2014. Muscle nutrients of the backcross progeny of female diploid F-1 hybrid (blunt snout bream ×
MA
topmouth culter) × male blunt snout bream and its parents. J. Fish. China 38, 1786-1792. Hu, J., Wen, H.S., Guan, J., Guan, S.G., He, F., Li, J.F., Shi, D., Ma, R.J., Liu, M., Mu, W.J., Zhang, Y.Q., 2012. Cloning of IGFBP-1, -2 and expression analys is during adult and early developmental stages in Scophthalmus maximus. Acta Oceanol. Sin. 34, 139-146.
ED
Hwa V, Oh Y, Rosenfeld RG (1999) The insulin-like growth factor-binding protein (IGFBP) superfamily. Endocr. Rev. 20, 761-787.
EP T
Jayaraman, G., Srimathi, S., Bjarnason, J.B., 2006. Conformation and stability of elastase from Atlantic cod, Gadus morhua. BBA-Gen Subjects. 1760, 47-54. Jia, Y. Y., Gu, Z.M., Ye, J.Y., Yang, Y.J., Zhu, J.J., Huang, X.M., 2011. Analys is on genetic variations of Erythroculter ilishaeformis (♀) × Megalobrama amblycephala (♂) hybrids F1 by SRAP markers. Journal of Shanghai Ocean University 20, 198-203.
AC C
Jiang, W.P., Jia, Y. Y., Liu, S.L., Li, Q., Li, T., Gu, Z.M., 2016. Comparative analysis of intermuscular bones in hybrid F1, F2 of (C. alburnus) (♀) × (M. amblycephala) (♂) and its parents. Acta Hydrobiologica Sinica 40, 277-286. Kalloohosein, H.E., Whitehead, J.P., Soos, M., Tavare, J.M., Siddle, K., Orahilly, S., 1997. Differential signaling to glycogen synthesis by the intracellular domain of the insulin versus the insulin-like growth factor-1 receptor. Evidence from studies of TrkC-chimeras. J. Biol. Chem. 272, 24325-24332. Kimball, S.R., Shantz, L.M., Horetsky, R.L., Jefferson, L.S., 1999. Leucine regulates translation of specific mRNAs in L6 myoblasts through mTOR-mediated changes in availability of eIF4E and phosphorylation of ribosomal protein S6. J. Biol. Chem. 274, 11647-11652. Koboldt, D.C., Chen, K., Wylie, T., Larson, D.E., McLellan, M.D., Mardis, E.R., Weinstock, G.M., Wilson, R.K., Ding, L., 2009. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics 25, 2283-2285.
21
ACCEPTED MANUSCRIPT Lafarga-De, La., Cruz, F., Nunez-Acuna, G., Gallardo-Escarate, C., 2013. Hybridization between Haliotis rufescens and Haliotis discus hannai: evaluation of fertilization, larval development, growth and thermal tolerance. Aquac. Res. 44, 1206-1220. Leng, X.J., Wu, X.F., Tian, J., Li, X.Q., Guan, L., Weng, D.C., 2012. Molecular cloning of fatty acid synthase from grass carp (Ctenopharyngodon idella) and the regulation of its expression by dietary fat level. Aquacult. Nutr. 18, 551-558. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. Bmc Bioinformatics 12, 323.
PT
Li, H., 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987-2993.
RI
Li, X.F., Jiang, Y.Y., Liu, W.B., Ge, X.P., 2012. Protein-sparing effect of dietary lipid in practical diets for blunt snout bream (Megalobrama amblycephala) fingerlings : effects on digestive
SC
and metabolic responses. Fish Physiol. Biochem. 38, 529-541.
Li, S.W., Wang, D., Cao, Y.S., Zhang, Y., Liu, H.B., Lu, T.Y., 2018. Transcriptome profile of amur sturgeon (Acipenser schrenckii) liver provides insights into immune modulation in response
NU
to Yersinia ruckeri infection. Aquaculture 492, 137-146.
Livak, K.J., Schmittgen, T.D., 2001. Analys is of relative gene expression data using real-time
MA
quantitative PCR and the 2(T)(-Delta Delta C) method. Methods 25, 402-408. Lu, Z.J., Xie, S.L., Wang, C., Zhou, A.G., Zou, J.X., 2015. Comparative study on the growing and digestion indexes of Snakehead with nine mating groups. Feed Industry 36, 25-30. Ma, X.M., Blenis, J., 2009. Molecular mechanisms of mTOR-mediated translational control. Nat.
ED
Rev. Mol. Cell Biol. 10, 307-318.
Mahardini, A., Yamauchi, C., Takeuchi, Y., Rizky, D., Takekata, H., Takemura, A., 2018. Changes
EP T
in mRNA abundance of insulin-like growth factors in the brain and liver of a tropical damselfish, Chrysiptera cyanea, in relation to seasonal and food-manipulated reproduction. Gen. Comp. Endocr. 269, 112-121. Morais, S., Monroig, O., Zheng, X., Leaver, M.J., Tocher, D.R., 2009. Highly unsaturated fatty acid synthesis in Atlantic salmon: characterization of ELOVL5- and ELOVL2-like elongases.
AC C
Mar. Biotechnol. 11, 627-639.
Moro, G.V., Camilo, R.Y., Moraes, G., Fracalossi, D.M., 2010. Dietary non-protein energy sources: growth, digestive enzyme activities and nutrient utilization by the catfish jundia, Rhamdia quelen. Aquac. Res. 41, 394-400. Ozsolak, F., Milos, P.M., 2011. RNA sequencing: advances, challenges and opportunities. Nat. Rev. Genet. 12, 87-98. Salze, G., Mclean, E., Craig, S.R., 2012. Dietary taurine enhances growth and digestive enzyme activities in larval cobia. Aquaculture 362-363, 44-49. Smith, S., Witkowski, A., Joshi, A.K., 2003. Structural and functional organization of the animal fatty acid synthase. Prog. Lipid Res. 42, 289-317. Stewart, C. E., Rotwein, P., 1996. Growth, differentiation, and survival: multiple physiological functions for insulin-like growth factors. Physiol. Rev. 76,1005-1026. Stitt, T.N., Drujan, D., Clarke, B. A., Panaro. F., Timofeyva, Y., Kline, W.O., Gonzalez, M.,
22
ACCEPTED MANUSCRIPT Yancopoulos, G.D., Glass, D.J., 2004. The IGF-1/PI3K/Akt pathway prevents expression of muscle atrophy-induced ubiquitin ligases by inhibiting FOXO transcription factors. Mol. Cell 14, 395-403. Sun, Y., Guo, C. Y., Wang, D. D., Li, X. F., Zhang, Y. H., You, X. X., Shi, Q., Hu, G. J., Fang, C., Lin, H.R., Zhang, Y., 2016a. Transcriptome analysis reveals the molecular mechanisms underlying growth superiority in a novel grouper hybrid (Epinephelus fuscogutatus♀ × E. lanceolatus♂). BMC Genetics 17, 24. Sun, Y., Huang, Y., Hu, G.J., Zhang, X.H., Ruan, Z.Q., Zhao, X.M., Guo, C.Y., Tang, Z.J., Li, X.F.,
PT
You, X.X., Lin, H.R., Zhang, Y., Shi, Q., 2016b. Comparative transcriptomic study of muscle provides new insights into the growth superiority of a novel grouper hybrid. Plos One 11, e0168802.
RI
Verma, E., Chakraborty, S., Tiwari, B., Mishra, A.K., 2018. Transcriptional regulation of acetyl CoA and lipid synthesis by PII protein in Synechococcus PCC 7942. J. Basic Microb. 58,
SC
187-197.
Wang, C.S., Hartsuck, J.A., 1993. Bile salt-activated lipase. A multiple function lipolyt ic enzyme. Biochim. Biophys. Acta. 1166, 1-19.
NU
Wang, S., Tang, C.C, Tao, M., Qin, Q.B., Zhang, C., Luo, K.K., Zhao, R.R., Wang, J., Ren, L., Xiao, J., Hu, F.Z., Zhou, R., Duan, W., Liu, S.J., 2019. Establishment and application of
MA
distant hybridization technology in fish. Sci. China Life Sci. 62, 22-45. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57-63.
Xiao, J., Kang, X.W., Xie, L.H., Qin, Q.B., He, Z.L., Hu, F.Z., Zhang, C., Zhao, R.R., Wang, J.,
ED
Luo, K.K., Liu, Y., Liu, S.J., 2014. The fertility of the hybrid lineage derived from female Megalobrama amblycephala × male Culter alburnus. Anim. Reprod. Sci. 151, 61-70.
EP T
Xie, C., Mao, X.Z., Huang, J.J., Ding, Y., Wu, J.M., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L.P., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316-W322. Xu, W.N., Qian, Y., Li, X.F., Li, J.Y., Li, P.F., Cai, D.S., Liu, W.B., 2017. Effects of dietary biotin on growth performance and fatty acids metabolism in blunt snout bream, Megalobrama
AC C
amblycephala fed with different lipid levels diets. Aquaculture 479, 790-797. Zhai, R.R., Feng, Y., Wang, H.M., Zhan, X.D., Shen, X.H., Wu, W.M., Zhang, Y.X., Chen, D.B., Dai, G.X., Yang, Z.L., Cao, L.Y., Cheng, S.H., 2013. Transcriptome analys is of rice root heterosis by RNA-Seq. Bmc Genomics 14, 19. Zhang, C., Lin, C., Fu, F., Zhong, X., Peng, B., Yan, H., Zhang, J., Zhang, W., Wang, P., Ding, X., Zhang, W., Ding, X., Zhang, W., Zhao, L., 2017. Comparative transcriptome analysis of flower heterosis in two soybean F1 hybrids by RNA-seq. PloS one 12, e0181061. Zhang, H. Y., He, H., Chen, L.B., Li, L., Liang, M.Z., Wang, X.F., Liu, X.G., He, G.M., Chen, R.S., Ma, L.G., Deng, X.W., 2008. A genome-wide transcription analysis reveals a close correlation of promoter INDEL polymorphism and heterotic gene expression in rice hybrids. Mol. Plant 1, 720-731. Zheng, G.D., Guo, D.D., Wu, C.B., Chen, J., Jiang, X.Y., Zou, S.M., 2019. The obvious heterosis and genetic characters of intergeneric cross and backcross juveniles between blunt snout
23
ACCEPTED MANUSCRIPT bream (Megalobrama amblycephala) and topmouth culter (Culter alburnus). Aquac. Res. 50, 1634-1643. Zheng, G.D., Wang, C.L., Guo, D.D., Jiang, X.Y., Zou, S.M., 2017. Ploidy level and performance in meiotic gynogenetic offsprings of grass carp using UV-irradiated blunt snout bream sperm. Aquaculture and Fisheries 2, 213-219. Zheng, G.D., Zhang, Q.Q., Li, F.G., Chen, J., Jiang, X.Y., Zou, S.M., 2015. Genetic characteris tics and growth performance of different Megalobrama amblycephala (♀) × Erythroculter
NU
SC
RI
PT
ilishaeformis (♂) hybrids. J. Fish. Sci. China 22, 402-409.
MA
Figure legends Fig. 1. Landscape of unigene distribution in the backcross BC. (A) E-value distribution of unigenes searched against public databases with an E-value cut-off of 1E-5. (B) Identity distribution of unigenes searched against public databases with an E-value cut-off of 1E-5. (C) Unigenes conserved in BC and model species. Unigenes of the BC were characterized by species by using BLASTX.
ED
searching against public databases. The number of BC homologous genes identified in other
EP T
Fig. 2. DEGs analysis and volcano plot for ‘BC vs. MA’ (A) and ‘BC vs. CA’ (B). The x-axis is the value of Log2 (Fold Change), and the y-axis is the value of Log2 (p-value). The red and yellow dots reveal the up-regulated DEGs, the blue and light blue dots reveal the down-regulated DEGs. MA, CA and BC denote M. amblycephala, C. alburnus and their backcrossing offsprings,
AC C
respectively.
Fig. 3. DEGs and the overlaps of DEGs between ‘BC vs. MA’ and ‘BC vs. CA’. Fig. 4. GO enrichment of DEGs in hepatopancreas of BC and their parents. (A) GO enrichment of DEGs of ‘BC vs. MA’. (B) GO enrichment of DEGs of ‘BC vs. CA’. Asterisk showed the most enriched GO terms in each category. Fig. 5. Hierarchical cluster analysis of DEGs involved in the growth superiority. The color key represents FPKM normalized log2 transformed counts in hepatopancreas of three fish species. Changes in expression levels are shown using color scales with saturation at > 2-fold changes. Red and green gradients indicate a increase and decrease of DEGs, respectively.
24
ACCEPTED MANUSCRIPT Fig. 6. A summary of the SSRs identified from the MA, CA and BC. Fig. 7. Distribution of putative single nucleotide polymorphisms (SNP) in the MA, CA and BC. Fig. 8. Validation of RNA-seq data of ten DEGs by qRT-PCR. β-actin was used as an internal control and each value represents average of three separate biological replicates. Fig. 9. The overview map of molecular mechanisms underlying growth superiority of the BC.
PT
Blue arrows denote the DEGs of ‘BC vs. MA’. Red arrows denote the DEGs ‘BC vs. CA’. Up or down arrows stand for up- or down-regulation in the BC compared with their parent. CM: cell membrane; MM: mitochondrial membrane; ERM: endoplasmic reticulum membrane; ACC: acetyl
RI
CoA carboxylase; oaa: oxaloacetate.
SC
Supplementary Information Supplementary Fig. 1. Length distribution of transcripts and unigenes in BC (A, D), MA (B, E)
NU
and CA (C, F) respectively.
Supplementary Fig. 2. COG functional classification of unigenes in BC. A total of 7,406
MA
assembled unigenes were annotated and assigned to 25 functional categories. Supplementary Fig. 3. KEGG pathway enrichment analysis of the DEGs of ‘BC vs. MA’. The x-axis is the significant enriched KEGG pathway classification and y-axis is P-value of
ED
enrichment.
enrichment.
EP T
Supplementary Fig. 4. KEGG pathway enrichment analysis of the DEGs of ‘BC vs. CA’. The x-axis is the significant enriched KEGG pathway classification and y-axis is P-value of
Supplementary Table 1. COG, KOG and NOG functional classification of unigenes in BC, MA
AC C
and CA, respectively.
25
ACCEPTED MANUSCRIPT Table 1. Genes and specific primers used for quantitative real-time PCR. Forward primer sequence (5’-3’)
Reverse primer sequence (5’-3’)
TRY
GGGTGTCTGATGACCTTAGTG
TGTCTGCTGCTCATTGCT
ELA1
TCTCCTCCTCCTTCAGGTTATG
CTTGTGGAGGCAGCCTTATT
CPA2
GCTTTCACTCACACCAATAACC
CCAGCATCCCAGTTCCTATT
BAL
GAGGAGATCGCAAAGAAGGTAG
ATCCACATTCCCAGCAAGAG
IGF1
CAAGAGAGGAGTTTGCAGTGA
GCAAGGGTTCCATCTGGTATAA
IGFBP2b
GCTGACAGAGGGCTAGATTTG
CTGCAGGCAAACTTGTGTTTAT
RAPTOR
CTGCCTCACTACACCAATCAA
GCCTGGGATCTTCTCAATCAA
EIF4E
CAGATGGGCTCTCTGGTATTTC
AGTTTACTGGGCTGCTGTATG
FAS
AGGTGTCCGAGAATGGAAATC
ELOVL1
CCGGCGTACAGTACAAAGAA
AGACCCAAGGTTGAAGGATTAC
β-actin
CGTGCTGTTTTCCCTTCCATT
CAATACCGTGCTCAAAGGATACTT
SC
RI
PT
Gene name
AC C
EP T
ED
MA
NU
GTCTGCAGGCATTGGTTTATTC
26
ACCEPTED MANUSCRIPT
Table 2. Growth performance and digestive enzyme activities (U g-1 tissue protein) Group
IW
FW
SGR
WGR
HSI
CF
Protease
Lipase
BC
6.50±0.35
50.70±3.47a
2.28a
680.00a
1.48±0.17ab
2.10±0.02a
31.60±1.34a
22.18±2.02a
MA
6.50±0.34
44.81±3.24b
2.14ab
589.38b
1.68±3.24a
1.86±0.01ab
23.50±1.10b
20.16±1.50b
CA
6.50±0.62
28.25±2.02c
1.63c
334.61c
1.07±0.11b
1.11±0.03c
25.43±1.15b
18.60±2.10b
PT
Notes: IW (g): initial body weight, FW (g): final body weight, Weight gain rate (WGR, %) = (FW − IW)×100 / IW, Specific growth rate (SGR, %) = (lnFW − lnIW)
RI
×100/number of culture days, Hepatosomatic index (HSI, %) = 100×hepatopancreas
SC
weight/body weight, Condition factor (CF, %) = 100×(body weight, g)/(body length, cm)3 . Digestive enzyme activities were expressed as micromoles of hydrolyzed
NU
substrate min-1 g-1 tissue protein (U g-1 tissue protein). Means in the same column with
AC C
EP T
ED
MA
different superscripts are significantly different (P < 0.05).
27
ACCEPTED MANUSCRIPT Table 3. Characterization of raw/clean data and mapping rate Total raw
# of clean
Total clean
Number of reads
Mapping
reads
read size (bp)
reads
read size (bp)
mapped to assembly
ratea (%)
BC
50,064,120
7,559,682,120
49,063,830
7,316,679,170
40,859,480
83.28
MA
56,209,354
8,487,612,454
55,225,738
8,244,390,148
49,405,514
89.46
CA
52,814,668
7,975,014,868
51,725,452
7,714,644,875
Average
53,029,381
8,007,436,481
52,005,007
7,758,571,398
Total
159,088,142
24,022,309,442
156,015,020 23,275,714,193
RI
SC
46,230,046
89.38
45,498,347
87.37
136,495,040
——
Reads with ≤ 5 base mismatches were counted when mapped to the reference
NU
a
PT
# of raw
EP T
ED
MA
sequences.
AC C
Sample
28
ACCEPTED MANUSCRIPT Table 4. Summary of cDNA sequences of the BC and their parents
MA
CA
BC
MA
CA
Total sequence No.
77102
63367
67087
62248
55182
57366
Total length (bp)
71136472
52919414
56482685
50596979
41110483
42828640
Largest length (bp)
16621
15752
15127
16621
15752
15127
Smallest length (bp)
201
201
201
201
201
201
Average length (bp)
923
835
842
813
745
747
N50
1799
1524
1601
1581
1288
1349
GC%
45.43
45.69
45.57
45.34
45.48
45.36
AC C
EP T
ED
MA
NU
SC
BC
PT
Unigenes
RI
Transcripts
29
ACCEPTED MANUSCRIPT
Gene name
BC-FPKM
MA-FPKM
CA-FKPM
Digestive
TRY
4614.88±65.81
411.07±6.83
1658.76±35.34
enzyme
ELA1
1327.78±25.27
145.72±5.47
400.73±11.26
CTRL1
564.14±18.29
6.51±0.84
125.65±11.31
CPA2
92.67±5.87
5.03±0.18
20.96±2.45
BAL
425.73±20.28
20.12±2.37
61.82±8.46
IGF1
28.02±1.05
2.73±0.84
3.11±0.93
IGF2a
82.29±5.86
27.62±3.63
13.79±2.86
IGFBP1
28.59±1.04
137.95±3.86
241.57±10.81
IGFBP2b
57.74±3.17
11.62±0.89
9.01±1.04
Protein
PI3KR
56.03±5.61
5.46±1.15
6.23±0.86
synthesis
RAPTOR
182.19±25.88
61.33±6.43
90.33±8.89
EIF4E
152.65±12.31
62.36±9.39
13.88±2.84
Fatty acid
CS
101.67±16.81
19.13±1.87
27.50±1.03
synthesis
MDH
979.76±32.12
404.63±15.23 469.19±21.13
1379.37±59.11
658.22±37.27 82.96±6.82
ELOVL1
33.86±5.45
14.14±4.72
6.65±1.37
ELOVL5
283.73±27.80
73.77±9.17
83.06±8.82
1213.94±55.88
604.20±30.28 35.71±3.91
AC C
EP T
FAS
ELOVL6
RI
SC
NU
ED
GH/IGF axis
PT
Function
MA
Table 5. FPKMs of digestive enzyme, IGF system and protein/fatty acid synthesis related genes in BC, MA and CA.
30
ACCEPTED MANUSCRIPT Highlights Intergeneric backcross BC with growth superiority and enhanced digestive enzyme activity were obtained between Megalobrama amblycephala and Culter alburnus.
PT
134 differentially expressed genes (DEGs) were identified between BC and
RI
their parents in hepatopancreas by transcriptome analysis.
SC
The DEGs related to the IGF system, digestive enzyme, protein/fatty acid
AC C
EP T
ED
MA
NU
synthesis may contribute to growth advantage of backcross BC.
31
Guodong Zheng, Chenbin Wu, Juan Liu, Jie Chen, Shuming Zou PII: DOI: Article Number: Reference:
S0044-8486(19)30779-3 https://doi.org/10.1016/j.aquaculture.2019.734317 734317 AQUA 734317
To appear in:
aquaculture
Received date: Revised date: Accepted date:
2 April 2019 30 June 2019 16 July 2019
Please cite this article as: G. Zheng, C. Wu, J. Liu, et al., Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂, aquaculture, https://doi.org/10.1016/j.aquaculture.2019.734317
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Transcriptome analysis provides new insights into the growth superiority of a novel backcross variety, Megalobrama amblycephala ♀ × (M. amblycephala ♀ × Culter alburnus ♂) ♂
PT
Guodong Zheng, Chenbin Wu, Juan Liu, Jie Chen, Shuming Zou*
RI
Key Laboratory of Freshwater Aquatic Genetic Resources, Ministry of Agriculture,
SC
Genetics and Breeding Center for Blunt Snout Bream, National Demonstration Center
NU
for Experimental Fisheries Science Education, Shanghai Ocean University, Shanghai 201306, China
MA
*Correspondence Authors: Dr. Shu-Ming Zou
Tel.: +86-021-61900345
ED
Fax: +86-021-61900345
AC C
EP T
Email: [email protected]
1
ACCEPTED MANUSCRIPT
Abstract A novel backcross variety (BC) [Megalobrama amblycephala (MA) ♀ × (Megalobrama amblycephala ♀ × Culter alburnus (CA) ♂) ♂] was generated by artificial insemination. The BC exhibits growth superiority and enhanced digestive
PT
enzyme activity than that of its parents. We conducted the transcriptome analysis
RI
using hepatopancreas between the BC and its parents. Three hundred and twenty-five
SC
and 872 differentially expressed genes (DEGs) were identified in ‘BC vs. MA’ and
NU
‘BC vs. CA’, respectively. The GO and KEGG enrichment analysis of DEGs were classified. From the 134 overlapping DEGs and their pathway, the digestive enzyme
MA
genes, such as TRY, ELA1, CTRL1, CPA2, and BAL, and the IGF system genes, such as IGF1, IGF2a, and IGFBP2b, of the BC had significant difference in expression.
ED
Furthermore, the expression of downstream genes related to the protein and fatty acid
EP T
synthesis pathways, such as PI3KR, RAPTOR, EIF4E, CS, MDH, FAS, ELOVL1, ELOVL5, and ELOVL6, were significantly up-regulated. In addition, a total of 40732
AC C
SSRs and 52951 SNPs were identified from the transcriptome. The results revealed that the enhanced digestive enzymes activity and synthesis of proteins/fatty acids might be contributed to growth superiority of backcross BC. Key words: transcriptome; backcross; growth superiority; digestive enzyme; protein/fatty acid synthesis
2
ACCEPTED MANUSCRIPT
1. Introduction Blunt snout bream (Megalobrama amblycephala, MA) and topmouth culter (Culter alburnus, CA) are freshwater fish species of the subfamily Cultrinae. In China, both species are important in the fish polyculture system. Although the M.
PT
amblycephala grows faster than the C. alburnus, but its meat quality is inferior to that
RI
of the C. alburnus. To produce new varieties with growth superiority and high quality
SC
meat, hybridization has been introduced into the artificial breeding of these two
NU
species. According to statistics, the yield of the hybrids between M. amblycephala and C. alburnus has exceeded 50 thousand tons in China in 2018 (China Aquaculture
MA
Network, www.shuichan.cc). In previous study, a fertlie hybrid lines between M. amblycephala and C. alburnus were established (Xiao et al., 2014), this hybrid were
ED
also found to have obvious heterosis in growth performance, feed utilization and
EP T
muscle quality (Zheng et al., 2015; Zheng et al., 2019). Moreover, the total amino acid and essential amino acid content of hybrids between M. amblycephala and C.
AC C
alburnus were higher than those of their parents (He et al., 2014), SRAP marker analysis showed that the genetic diversity of this hybrids increased significantly (Jia et al., 2011) and the number of intermuscular bones decreased significantly (Jiang et al., 2016). In addition, a similar hybrid between M. terminalis and C. alburnus also showed growth advantage (Guo et al., 2018). RNA-Seq has been widely used as an effective tool for transcriptome analysis in order to discover, profile, and quantify the RNA transcripts (Wang et al., 2009).
3
ACCEPTED MANUSCRIPT
RNA-Seq combines the advantages of high throughput, low background noise and high sensitivity, making it feasible to map transcribed regions, quantify gene expression, and distinguish different isoforms and allelic expression (Wang et al., 2009; Ozsolak and Milos, 2011). Recently, the differently expressed genes (DEGs)
PT
and pathways related to heterosis have been identified through the transcriptome
RI
profiles analysis in some commercially crops and fishes, such as hybrid maize (Zea
SC
mays L.) (Bi et al., 2014), hybrid rice (Oryza sativa) (Zhai et al., 2013), hybrid
NU
soybean (Glycine max (L.) Merr) (Zhang et al., 2017), hybrid pufferfish (Takifugu) (Gao et al., 2013), hybrid grouper (Epinephelus spp.) (Sun et al., 2016a, b), etc. These
MA
studies could help us to further understand the mechanisms underlying the heterosis in crops and fishes.
ED
In order to reveal the molecular mechanism of the heterosis of M. amblycephala
EP T
and C. alburnus, transcriptome analysis was carried out for the two species and their hybrids. It is hoped that effective DEGs can be screened to reasonably explain the
AC C
heterosis of fish, and a large number of molecular markers can be developed for the study of QTL analysis and genetic maps.
2. Materials and methods 2.1. Experimental fish and feeding trial M. amblycephala (MA) and C. alburnus (CA) broodstocks were obtained from the Genetics and Breeding Center for Blunt Snout Bream of Ministry of Agriculture,
4
ACCEPTED MANUSCRIPT Shanghai, China. MC hybrids (MA♀ × CA♂) were produced first, when MC hybrids reached sexual maturity, the backcross BC (MA♀ × MC♂) was produced and parental selfbred strains of MA (♀ × ♂) and CA (♀ × ♂) were also produced by artificial insemination. All fish used in the study came from a single cross family, eggs from
PT
the same female used for the backcross with the male hybrid also used for maintaining
RI
the parental strains. These three crosses (BC, MA and CA) with the same age and
SC
similar sizes were randomly distributed into 3 concrete ponds (6 × 4 × 1.2 m, L : W :
NU
H) with water depth of 1m. Each pond was stocked all three crosses and each cross had 50 fry (total of 150 fish per pond). The feeding period was 90 days. During the
MA
feeding trial, fish were fed twice (06:30 and 17:30) each day with equal amounts of floating compound feed (containing 33% protein, 3% lipid and 8.5% fiber, Tech-bank
ED
Co., Ltd, Ningbo, China) and until the feed is eaten completely by the fish. One half
EP T
of the pond water was changed weekly to ensure fresh water quality. Temperature of
AC C
water ranged from 22~30 ℃ and dissolved oxygen was above 5.0 mg/L.
2.2. Growth measurement and sample collection At the end of feeding trial, after 24 hours of starvation, the weight gain rate (WGR), specific growth rate (SGR), hepatosomatic index (HSI) and condition factor (CF) were calculated. Ten individuals per cross were euthanized with 100 mg/L of MS222 (Sigma, USA) before tissues were collected. Hepatopancreases (for RNA extraction, sequencing and digestive enzyme activity assays) was collected and stored
5
ACCEPTED MANUSCRIPT
immediately at -80°C.
2.3. Digestive enzyme activity assays The whole hepatopancreas of three samples per cross was weighed and
PT
homogenized (dilution 1:10) in ice-cold buffer. The homogenization buffer solution
RI
was 0.02 M Tris/0.01 M phosphate, pH 7.0, in 50% (v/v) glycerol (Moro et al., 2010).
SC
The extract was centrifuged at 3000 rpm for 10 min at 4°C, and the supernatant was
NU
used as the enzyme source using BSA (Sigma, USA) as standard. The activity of total protease was detected using Folin-phenol method and the activity of total lipase was
MA
detected using commercial kits (ref. no. A054) produced by Jiancheng Bioengineering
ED
Institute (Nanjing, China) (Li et al., 2012).
EP T
2.4. RNA extraction, library construction and sequencing Total RNA from hepatopancreas samples was extracted using TRIzol reagent
AC C
(Invitrogen, Carlsbad, CA, USA) and purified using the RNeasy Mini Kit (Qiagen, Valencia, CA, USA). RNA degradation and contamination were detected by 1% agarose gels. The integrity and concentration of RNA samples were examined by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Equal amounts of nine RNA samples (three samples per cross) were used for cDNA synthesis and RNA-seq. Before the construction of library, Magnetic Oligo-dT beads (Invitrogen, USA) were used to isolate poly (A) + mRNA, and then cDNA
6
ACCEPTED MANUSCRIPT libraries were constructed and sequenced by Illumina Hiseq™ 4000 platform (Majorbio Biotech Co., Ltd, Shanghai, China). Briefly, 8 μg of total RNA for each sample was used to construct libraries by using Truseq TM RNA sample prep Kit (Illumina, USA) according to the manufacturer’s protocol. The constructed DNA
PT
template was enriched by PCR amplification (15 cycles) . Amplicons were collected
RI
and purified by Certified Low Range Ultra Agarose (Bio-Rad, USA) gel
SC
electrophoresis. Before sequencing, the DNA libraries were quantified by using
NU
TBS-380 micro fluorometer with Picogreen® reagent (Invitrogen, USA). Clone clusters were generated on Illumina cBot, using Truseq PE Cluster Kit v3-cBot-HS,
MA
and high- throughput sequencing was performed on an Illumina Hiseq™ 4000
ED
sequencer, using Truseq SBS Kit (300 cycles).
EP T
2.5. De novo assembly and functional annotation of unigenes The raw paired end reads from all transcriptomes were cleaned by removing
AC C
adapter contamination, low quality sequences (reads with over 10% unknown base pairs‘N’) and empty reads. All clean reads were then assembled using the denovo assembly program Trinity. The assembled unigenes of three crosses were used for BlastX search and annotation against the NR, Swissprot, String (Franceschini et al., 2013), COG (Clusters of Orthologous Groups of proteins), KEGG (Kyoto Encyclopedia of Genes and Genomes) databases with an E-value cut-off of 10−5 . Functional annotation including biological process, molecular function and cellular
7
ACCEPTED MANUSCRIPT
component for gene ontology (GO) terms was analysed by using BLAST2GO software (Conesa et al., 2005).
2.6. Comparative expression analysis
PT
The read counts were further normalized into FPKM (Fragments Per Kilobase of
RI
exon model per Million mapped reads) values. The FPKM values from the three
SC
libraries were pairwise compared, and the fold changes (FC) of the genes were
NU
calculated by using RSEM software (version 1.2.7) (Li and Dewey, 2011). If |log2(FC)| ≥ 1 and false discovery rate (FDR) ≤ 0.05, these genes were considered as
MA
having significant differential expression (Anders and Huber, 2010). The expressions of unigenes were displayed intuitively by Volcano Plot. After DEGs in the
ED
hepatopancreas between backcrosses and their parents were screened, the enrichment
EP T
analysis of GO and KEGG pathways was carried out by Goatools software (version0.4.7) and KOBAS 2.0 software (Xie et al., 2011). The R program was used
AC C
to depict the heat maps for gene clustering in the present study.
2.7. Detection of SSRs and SNPs Assembled
transcriptomes
were
screened
for
microsatellites
using
Msatcommander program (Faircloth, 2008). The parameters were designed to identify perfect mono-hexa- nucleotide motifs with a minimum of ten repetitions for mononucleotides, six repetitions for dinucleotides, and five repetitions for tri- hexa
8
ACCEPTED MANUSCRIPT nucleotides. Potential SNPs were filtered using the program Samtools (Li, 2011) and VarScan v.2.2.7 (Koboldt et al., 2009 ) with the default parameters.
2.8. Qualitative real-time PCR (qRT-PCR)
PT
Total RNA of hepatopancreas was extracted with the method described above, and
RI
cDNA was synthesized by PrimeScript RT reagent Kit (Takara, Janpan). The cDNA
SC
was then used for qRT-PCR analysis using specific primers (Table 1). qRT-PCR was
NU
performed by using SYBR Green RT-PCR kit (Takara, Janpan) on the CFX96 Touch™ real- time PCR Detection System (BioRad, USA). The housekeeping gene
MA
β-actin had been proved to be stable between the comparison groups and then used as the internal control. All reactions were performed in triplicates. The relative
ED
expression ratio of the target genes versus β-actin gene was calculated using 2-ΔΔCT
EP T
method (Livak and Schmittgen, 2001).
AC C
2.9 Statistical analysis
Digestive enzyme activity and gene expression data were carried out by using SPSS version 17 (Michigan Avenue, Chicago, IL, USA) for significant differences. If significant (P < 0.05) differences were found in factors, Duncan's multiple range test (Duncan, 1955) was used to rank the means.
3. Result
9
ACCEPTED MANUSCRIPT
3.1. Growth performance and digestive enzyme activity The growth performances of the backcross BC, MA, and CA are presented in Table 2. The weight gain rate (WGR) of BC, MA, and CA was 680.00%, 589.38%, and 334.61%, respectively. The WGR of BC was 15.38% and 103.22% higher than
PT
that of its parents MA and CA, this result showed that the growth superiority of the
RI
BC was significant (P < 0.05). The condition factor (CF) of BC was slightly and
SC
significantly (P < 0.05) higher than that of MA and CA, respectively. However,
NU
hepatosomatic index (HSI) of BC was intermediated between their two parents. In addition, the protease and lipase activity in hepatopancreas of BC was significantly (P
MA
< 0.05) higher than that of their parents, protease activity of CA was slightly higher
ED
than that of MA and lipase activity of MA was higher than that of CA (P < 0.05).
EP T
3.2. High-throughput sequencing and mapping to the reference genome Over 7.45 Gb of raw data for each sample and 159 million 100-bp paired-end
AC C
reads with an average of 53 million reads for each of the three samples (Table 3) were obtained. The raw data were cleaned and quality checked, and then assembled. Approximately, 49 million (BC), 55 million (MA), and 51 million (CA) clean reads were produced; 83.28% (BC), 89.64% (MA), and 89.38% (CA) of these clean reads were mapped to the assembled sequences (≤ 5 base mismatches) (Table 3). An assembly of the reads generated 77102 (BC), 63367 (MA) and 67087 (CA) transcripts with an average length of 923 bp, 835 bp and 842 bp, respectively. And a
10
ACCEPTED MANUSCRIPT
total of 62248 (BC), 55182 (MA) and 57366 (CA) unigenes were generated, with an average length of 813 bp, 745 bp and 747 bp, respectively (Table 4). An average of 61.78% of the transcripts and 66.29% of the unigenes were < 600 bp in size for the three samples (Supplementary Fig. 1). The results showed that ORFs of 41713
PT
unigenes (67.01%) in BC, 34654 unigenes (62.79%) in MA and 34602 unigenes
RI
(60.32%) in CA were successfully predicted using Trinity respectively. Overall, there
SC
were 25 classifications produced from the COG/KOG/NOG-annotated putative
NU
proteins, including ‘metabolism’, ‘cellular processes a nd signaling’ and ‘information storage and processing’, in accordance with the GO annotation categories
MA
(Supplementary Fig. 2 and Supplementary Table 1).
Characterization of BC unigenes by searching against public databases: E-value
ED
distribution of the top hits in the databases showed 11205 (35.39%) unigenes with a
EP T
strong homology (<1.0e-50) (Fig. 1A), 27629 (87.28%) unigenes had a similarity higher than 80%, while 36754 (11.61%) showed a similarity between 60% and 80%
AC C
with respect to the identity distribution pattern. Therefore, 98.89% of the unigenes showing an identity higher than 60% along with a high-quality similar distribution supported the reliability of the de novo assembly (Fig. 1B). A total of 31657 putative known unigenes were found in all reference species, 16903 (53.39%) were found in Danio rerio (Fig. 1C).
3.3. Identification and analysis of the diff erently expressed genes (DEGs)
11
ACCEPTED MANUSCRIPT
FPKMs of each gene in the hepatopancreas of the backcross BC were compared with those of their parents. The DEGs of ‘BC vs. MA’ and ‘BC vs. CA’ were identified by filtering the data with the following criteria: absolute log2 (FC) ≥ 1, FDR < 0.05, and FPKM > 2. According to these criteria, volcano plots were used to
PT
describe the number of significant up and downregulated DEGs of ‘BC vs. MA’ (Fig.
RI
2A) and ‘BC vs. CA’ (Fig. 2B).
SC
As shown in Fig. 3, Venn diagram showed that there were 325 DEGs (130
NU
up-regulated and 195 down-regulated) in ‘BC vs. MA’ and 872 DEGs in (388 up-regulated and 484 down-regulated) ‘BC vs. CA’. Among those genes, 134
MA
overlapping DEGs were found in ‘BC vs. MA’ and ‘BC vs. CA’. Based on GO annotation, the hepatopancreas DEGs of ‘BC vs. MA’ and ‘BC vs. CA’ were enriched
ED
into 46 and 53 GO terms respectively, which provides an overview of ontology
EP T
content. Metabolic process, cell part and binding were the most three enriched GO terms in the biological process, cellular component and molecular function categories,
AC C
respectively (Fig. 4A and B). By KEGG analysis, the DEGs of ‘BC vs. MA’ were classified to 5 categories with 54 pathway and DEGs of ‘BC vs. CA’ were classified to 5 categories with 51 pathway (Supplementary Fig. 3-4), in which the protein digestion-absorption pathway and pancreatic secretion pathway were the most significant pathway in both ‘BC vs. MA’ and ‘BC vs. CA’.
3.4. Analysis of diff erently expressed genes related to growth
12
ACCEPTED MANUSCRIPT
By analyzing these 134 overlapping DEGs and their pathway, we found that 5 genes, including trypsinogen (TRY), pancreatic elastase 1 (ELA1), chymotrypsin- like 1 (CTRL1), carboxypeptidase 2 (CPA2) and bile salt-activated lipase (BAL) of BC had higher FPKMs value than both MA and CA, which provides further evidence that
PT
digestive enzyme genes were more active in the BC (Fig. 5 and Table 5).
RI
From these 134 overlapping DEGs and their pathway, we also analyzed the
SC
insulin- like growth factor (IGF) system related genes in hepatopancreases, such as
NU
IGF1, IGF2a, IGF binding protein1 (IGFBP1) and IGFBP2b were differentially expressed in BC when compared with those of its parents (Fig. 5 and Table 5). In
including
MA
addition, the DEGs were identified in protein and fatty acid synthesis pathways, phosphoinositide-3-kinase
regulatory
subunit
(PI3KR),
regulatory
ED
associated protein of mTOR (RAPTOR), eukaryotic translation initiation factor 4E
EP T
(EIF4E), citrate synthase (CS), malate dehydrogenase (MDH), fatty acid synthetase (FAS) and elongation of very long chain fatty acids protein (ELOVL1, ELOVL5 and
AC C
ELOVL6) (Fig. 5 and Table 5).
3.5. SSRs and SNP discovery Using the Msatcommander program, a total of 12718, 12922 and 15092 potential SSRs contained repeats of more than one nucleotides with a minimum of five repetitions were identified from MA, CA and BC, respectively (Fig. 6). The number of SSRs with ten repetitions is the largest in all three crosses and the BC group has
13
ACCEPTED MANUSCRIPT
more SSRs than MA and CA group. Among these SSRs, GT/GT, GAT and GATA were most common in dinucleotide, trinucleotide and quadnucleotide SSRs, respectively. A total of 52951 candidate SNPs were identified in 3 crosses. Of these SNP
PT
candidates, 38787 SNPs were putative transitions (Ts), and 14164 SNPs were putative
RI
transversions (Tv) with a mean Ts: Tv ratio of 2.73 (Fig. 7). The SNPs were then
SC
categorized into 4 classes, including class 1 (C/A, A/C, T/G and G/T) for 11.44%,
NU
class 2 (C/T, G/A, T/C and A/G) for 73.25%, class3 (C/G and G/C) for 6.84%, and
MA
class 4 (A/T and T/A) for 8.47%.
3.6. Validation of the RNA-seq analysis by qRT-PCR
ED
To validate the DEGs in FPKM values, ten DEGs were randomly selected for
EP T
qRT-PCR assay in hepatopancreas of these three crosses. The results were further compared with the FPKM values generated from RNA-seq. Our results showed that
AC C
the data of qRT-PCR were consistent with those of RNA-seq (Fig. 8).
4. Discussion
Hybridization is widely used to improve growth rate or other economic traits of animals and plants (Lafarga-De la Cruz et al., 2013; Zheng et al., 2017; Chen et al., 2018; Wang et al., 2019; Zheng et al., 2019). In recent years, many studies have focused on exploring the molecular basis of heterosis in hybrid species by
14
ACCEPTED MANUSCRIPT
transcriptome analysis (Zhang et al., 2008; Gao et al., 2013; Zhai et al., 2013; Sun et al., 2016a, b). In the present study, comparative transcriptomic analysis was performed on hepatopancreas of the backcross BC and their parents MA and CA. Among the annotated transcripts, 1197 DEGs (325 and 872 DEGs in ‘BC vs. MA’ and
PT
‘BC vs. CA’, respectively) were identified, and 134 overlapping DEGs were
RI
confirmed. By analyzing these 134 overlapping DEGs and their pathway, many
SC
growth superiority-related genes were found to be diff erentially expressed in BC
NU
compared with MA and CA. Based on the function of these DEGs and their location in the pathway, a overview map of the generation of growth superiority was produced
MA
(Fig. 9).
ED
4.1. Gene expression and activity of digestive enzyme
EP T
The growth and development of fish depends on nutrients, especially amino acids and fatty acids. High-quality digestive juices can break down food efficiently,
AC C
improving the feed utilization of fish (Wang and Hartsuck, 1993; Jayaraman et al., 2006). In the present study, we found that the expression of four protein digestive enzyme genes, such as TRY, ELA1, CTRL1 and CPA2, and one fat digestive enzyme gene BAL was higher in the BC than that of their parents (Fig. 9). Furthermore, we found that the activity of most protease and lipase also increased significantly in BC, which was consistent with digestive enzyme genes expression above. Some studies showed that the growth performance of fish was related to digestive enzyme activity
15
ACCEPTED MANUSCRIPT
(Salze et al., 2012; Lu et al., 2015). With increased activity of digestive enzymes, large molecules of protein and fat can be broken down more efficiently into small molecules of amino acids and fatty acids, which can provide more raw materials for
PT
protein/fatty acid synthesis.
RI
4.2. IGF system and protein synthesis
SC
The growth performance of vertebrates or fish is primarily effected by the
NU
growing axis (Duan 1997; Duan 1998). In this axis, IGF system is involved in cell regeneration, proliferation and differentiation (Stewart and Rotwein 1996; Mahardini
MA
et al., 2018). In the present study, four DEGs associated with IGF system were found between BC and their parents (Fig. 9). Among them, the expression of IGF1 and
ED
IGF2a in the hepatopancreas of BC was significantly up-regulated than that of their
EP T
parents, their main function was to amplify and activate downstream pathways (such as protein synthesis pathway mentioned in next paragraph) (Sun et al., 2016a).
AC C
However, IGFBP1 exhibited a significantly decreased expression in the BC compared with their parents. This might indicate that reduced IGFBP1 protein allow IGFs to be released from the IGFBP-IGFs conjugates, therefore, more separated IGFs could combine with IGFR on cell surface and exert their role, which is similar to previous study in Scophthalmus maximu (Hu et al., 2012). On the contrary, the expression of the IGFBP2b mRNA in the BC was up-regulated significantly, the same results were also found in hybrid grouper (Sun et al., 2016a).
16
ACCEPTED MANUSCRIPT
In present study, besides the increased expression of IGF1 and IGF2a, the genes PI3KR, RAPTOR and EIF4E in the protein synthesis pathway also had higher expression in the BC, indicating the strengthened ability of protein synthesis in the BC. This results are consistent with previous study of hybrid grouper (Sun et al.,
PT
2016a). Combined with the above paragraph, many phosphorylation events, including
RI
PI3K/AKT pathway, are activated by the binding of IGF-1 to its receptor IGFR1 (Stitt
SC
et al., 2004), thus leading to a series of anabolic effects, such as protein synthesis
NU
(Burgos and Cant, 2010) and glycogen synthesis (KallooHosein et al., 1997). When stimulated by IGF-1, the PI3K/AKT signal activates mTOR (Bibollet-Bahena and
MA
Almazan, 2009), mTORC1 is then formed by the combination of Raptor with mTOR, which can phosphorylate EIF4E-binding protein 1 (EIF4E-BP1) (Ma and Blenis,
ED
2009). The phosphorylated EIF4E-BP1 can promote the dissociation of EIF4E from
EP T
the inhibitory complex of EIF4E and EIF4E-BP1 (Kimball et al., 1999). Finally, the
AC C
increased available EIF4E improves the synthesis of protein (Kimball et al., 1999).
4.3. Fatty acid synthesis The expression of many genes involved in the fatty acid synthesis pathway, such as CS, MDH, FAS, ELOVL1, ELOVL5 and ELOVL6 was up-regulated in the hepatopancreas of backcross BC compared with those in their parents (Fig. 9). These up-regulated genes may lead to faster rate of fatty acid synthesis and their specific roles in fatty acid synthesis pathway were explained as follows. As the direct raw
17
ACCEPTED MANUSCRIPT
material for fatty acid synthesis, acetyl CoA can bind with oxaloacetate to generate citric acid by the catalysis of citrate synthetase (CS), and enter the cytoplasm from mitochondrion (Verma et al., 2018; Xu et al., 2017). Citric acid will be cleaved into acetyl CoA and oxaloacetate in cytoplasm. Oxaloacetate then has to enter the
PT
mitochondrion by the catalysis of MDH for the next cycle of combining acetyl CoA.
RI
In cytoplasm, acetyl CoA can turn into malonyl Co A by catalysis of acetyl CoA
SC
carboxylase (ACC), fatty acid synthetase (FAS) then catalyzes the transformation of
NU
both acetyl CoA and malonyl CoA into fatty acid (16C) (Smith et al., 2003; Leng et al., 2012). Additionally, the ELOVLs conduct the extension of very long chain fatty
MA
acid (18C or 20C or 24C) in endoplasmic reticulum (ER) to play a variety of
EP T
4.4. SSRs and SNPs
ED
biological roles (Morais et al., 2009; Gregory et al., 2010; Green and Olson, 2011).
By transcriptome analysis, a large number of SSRs and SNPs have been identified
AC C
in BC and their parents. We found that SSRs with 10 repeats and SNPs of transition (Ts) were the most numerous, respectively (Li et al., 2018). However, the number of SSRs and SNPs in BC was significantly higher than that of their parents, this might be because the backcross BC genome fused the genomes of both parents. This SSRs and SNPs will provide a helpful resource for marker-assisted selection (MAS), quantitative trait locus (QTL) association and population genetics analysis (Li et al., 2018).
18
ACCEPTED MANUSCRIPT
5. Conclusion In the present study, we performed a hepatopancreas transcriptome analysis for a backcross BC and their parents to explore the molecular mechanisms of growth
PT
superiority. The results revealed that the digestive enzymes activity and the synthesis
RI
of proteins/fatty acids might be important factors for growth superiority of BC.
SC
Moreover, a large number of SSRs were identified in hepatopancreas transcriptome,
NU
which would be beneficial to QTL analysis and genetic linkage. Overall, these
Acknowledgments
ED
growth superiority of hybrid fish.
MA
findings provide new insights into molecular and physiological basis underlying the
EP T
This work was supported by grants from the National Natural Sc ience Foundation of China (31572220), Project funded by China Postdoctoral Science Foundation
AC C
(2019M651473), and the Shanghai University Knowledge Service Platform (ZF1206).
19
SC
RI
PT
ACCEPTED MANUSCRIPT
NU
References Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome
MA
Biol. 11, R106.
Bi, Y.M., Meyer, A., Downs, G.S., Shi, X.J., El-Kereamy, A., Lukens, L., Rothstein, S.J., 2014. High throughput RNA sequencing of a hybrid maize and its parents shows different mechanisms responsive to nitrogen limitation. Bmc Genomics 15, 77.
ED
Bibollet-Bahena, O., Almazan, G., 2009. IGF-1-stimulated protein synthesis in oligodendrocyte progenitors requires PI3K/mTOR/Akt and MEK/ERK pathways. J. Neurochem. 109,
EP T
1440-1451. Burgos, S.A., Cant, J.P., 2010. IGF-1 stimulates protein synthesis by enhanced signaling through mTORC1 in bovine mammary epithelial cells. Domest. Anim. Endocrin. 38, 211-221. Chen, J., Luo, M., Li, S.N., Tao, M., Ye, X.L., Duan, W., Zhang, C., Qin, Q.B., Xiao, J., Liu, S.J., 2018. A comparative study of distant hybridization in plants and animals. Sci. China Life Sci.
AC C
61, 285-309.
Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analys is in functional genomics research. Bioinformatics 21, 3674-3676. Duan, C.M., 1997. The insulin-like growth factor system and its biological actions in fish. Am. Zool. 37, 491-503. Duan, C.M., 1998. Nutritional and developmental regulation of insulin-like growth factors in fish. J. Nutr .128, 306s-314s. Faircloth, B.C., 2008. Msatcommander: detection of microsatellite repeat arrays and automated, locus-specific primer design. Mol. Ecol. Resour. 8, 92-94. Firth, S.M., Baxter, R.C., 2002. Cellular actions of the insulin-like growth factor binding proteins. Endocr. Rev. 23, 824-854. Franceschini, A., Szklarczyk, D., Frankild, S., Kuhn, M., Simonovic, M., Roth, A., Lin, J.Y.,
20
ACCEPTED MANUSCRIPT Minguez, P., Bork, P., Von Mering C., Jensen, L.J., 2013. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 41, D808-D815. Gao, Y., Zhang, H., Gao, Q., Wang, L., Zhang, F., Siva, V.S., Zhou, Z., Song, L., Zhang, S., 2013. Transcriptome analysis of artificial hybrid pufferfish Jiyan-1 and its parental species: implications for pufferfish heterosis. PLoS One 8, e58453. Green, C.D., Olson, L.K., 2011. Modulation of palmitate-induced endoplasmic reticulum stress and apoptosis in pancreatic beta-cells by stearoyl-CoA desaturase and Elovl6. Am. J. Physiol.
PT
Endocrinol. Metab. 300, E640-649. Gregory, M.K., See, V.H., Gibson, R.A., Schuller, K.A., 2010. Cloning and functional characterisation of a fatty acyl elongase from southern bluefin tuna (Thunnus maccoyii).
RI
Comp. Biochem. Phys. B 155, 178-185.
Guo, H.H., Zheng, G.D., Wu, C.B., Chen, J., Jiang, X.Y., Zou, S.M., 2018. Comparative analys is
SC
of the growth performance and intermuscular bone traits in F 1 hybrids of black bream (Megalobrama terminalis) (♀) × topmouth culter (Culter alburnus) (♂). Aquaculture 492, 15-23.
NU
He, Z.L., Liu, S.J., Xiao J., Hu, F.Z., Wen, M., Ye, L.H., Zhang, C., Xu, K., Tao, M., 2014. Muscle nutrients of the backcross progeny of female diploid F-1 hybrid (blunt snout bream ×
MA
topmouth culter) × male blunt snout bream and its parents. J. Fish. China 38, 1786-1792. Hu, J., Wen, H.S., Guan, J., Guan, S.G., He, F., Li, J.F., Shi, D., Ma, R.J., Liu, M., Mu, W.J., Zhang, Y.Q., 2012. Cloning of IGFBP-1, -2 and expression analys is during adult and early developmental stages in Scophthalmus maximus. Acta Oceanol. Sin. 34, 139-146.
ED
Hwa V, Oh Y, Rosenfeld RG (1999) The insulin-like growth factor-binding protein (IGFBP) superfamily. Endocr. Rev. 20, 761-787.
EP T
Jayaraman, G., Srimathi, S., Bjarnason, J.B., 2006. Conformation and stability of elastase from Atlantic cod, Gadus morhua. BBA-Gen Subjects. 1760, 47-54. Jia, Y. Y., Gu, Z.M., Ye, J.Y., Yang, Y.J., Zhu, J.J., Huang, X.M., 2011. Analys is on genetic variations of Erythroculter ilishaeformis (♀) × Megalobrama amblycephala (♂) hybrids F1 by SRAP markers. Journal of Shanghai Ocean University 20, 198-203.
AC C
Jiang, W.P., Jia, Y. Y., Liu, S.L., Li, Q., Li, T., Gu, Z.M., 2016. Comparative analysis of intermuscular bones in hybrid F1, F2 of (C. alburnus) (♀) × (M. amblycephala) (♂) and its parents. Acta Hydrobiologica Sinica 40, 277-286. Kalloohosein, H.E., Whitehead, J.P., Soos, M., Tavare, J.M., Siddle, K., Orahilly, S., 1997. Differential signaling to glycogen synthesis by the intracellular domain of the insulin versus the insulin-like growth factor-1 receptor. Evidence from studies of TrkC-chimeras. J. Biol. Chem. 272, 24325-24332. Kimball, S.R., Shantz, L.M., Horetsky, R.L., Jefferson, L.S., 1999. Leucine regulates translation of specific mRNAs in L6 myoblasts through mTOR-mediated changes in availability of eIF4E and phosphorylation of ribosomal protein S6. J. Biol. Chem. 274, 11647-11652. Koboldt, D.C., Chen, K., Wylie, T., Larson, D.E., McLellan, M.D., Mardis, E.R., Weinstock, G.M., Wilson, R.K., Ding, L., 2009. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics 25, 2283-2285.
21
ACCEPTED MANUSCRIPT Lafarga-De, La., Cruz, F., Nunez-Acuna, G., Gallardo-Escarate, C., 2013. Hybridization between Haliotis rufescens and Haliotis discus hannai: evaluation of fertilization, larval development, growth and thermal tolerance. Aquac. Res. 44, 1206-1220. Leng, X.J., Wu, X.F., Tian, J., Li, X.Q., Guan, L., Weng, D.C., 2012. Molecular cloning of fatty acid synthase from grass carp (Ctenopharyngodon idella) and the regulation of its expression by dietary fat level. Aquacult. Nutr. 18, 551-558. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. Bmc Bioinformatics 12, 323.
PT
Li, H., 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987-2993.
RI
Li, X.F., Jiang, Y.Y., Liu, W.B., Ge, X.P., 2012. Protein-sparing effect of dietary lipid in practical diets for blunt snout bream (Megalobrama amblycephala) fingerlings : effects on digestive
SC
and metabolic responses. Fish Physiol. Biochem. 38, 529-541.
Li, S.W., Wang, D., Cao, Y.S., Zhang, Y., Liu, H.B., Lu, T.Y., 2018. Transcriptome profile of amur sturgeon (Acipenser schrenckii) liver provides insights into immune modulation in response
NU
to Yersinia ruckeri infection. Aquaculture 492, 137-146.
Livak, K.J., Schmittgen, T.D., 2001. Analys is of relative gene expression data using real-time
MA
quantitative PCR and the 2(T)(-Delta Delta C) method. Methods 25, 402-408. Lu, Z.J., Xie, S.L., Wang, C., Zhou, A.G., Zou, J.X., 2015. Comparative study on the growing and digestion indexes of Snakehead with nine mating groups. Feed Industry 36, 25-30. Ma, X.M., Blenis, J., 2009. Molecular mechanisms of mTOR-mediated translational control. Nat.
ED
Rev. Mol. Cell Biol. 10, 307-318.
Mahardini, A., Yamauchi, C., Takeuchi, Y., Rizky, D., Takekata, H., Takemura, A., 2018. Changes
EP T
in mRNA abundance of insulin-like growth factors in the brain and liver of a tropical damselfish, Chrysiptera cyanea, in relation to seasonal and food-manipulated reproduction. Gen. Comp. Endocr. 269, 112-121. Morais, S., Monroig, O., Zheng, X., Leaver, M.J., Tocher, D.R., 2009. Highly unsaturated fatty acid synthesis in Atlantic salmon: characterization of ELOVL5- and ELOVL2-like elongases.
AC C
Mar. Biotechnol. 11, 627-639.
Moro, G.V., Camilo, R.Y., Moraes, G., Fracalossi, D.M., 2010. Dietary non-protein energy sources: growth, digestive enzyme activities and nutrient utilization by the catfish jundia, Rhamdia quelen. Aquac. Res. 41, 394-400. Ozsolak, F., Milos, P.M., 2011. RNA sequencing: advances, challenges and opportunities. Nat. Rev. Genet. 12, 87-98. Salze, G., Mclean, E., Craig, S.R., 2012. Dietary taurine enhances growth and digestive enzyme activities in larval cobia. Aquaculture 362-363, 44-49. Smith, S., Witkowski, A., Joshi, A.K., 2003. Structural and functional organization of the animal fatty acid synthase. Prog. Lipid Res. 42, 289-317. Stewart, C. E., Rotwein, P., 1996. Growth, differentiation, and survival: multiple physiological functions for insulin-like growth factors. Physiol. Rev. 76,1005-1026. Stitt, T.N., Drujan, D., Clarke, B. A., Panaro. F., Timofeyva, Y., Kline, W.O., Gonzalez, M.,
22
ACCEPTED MANUSCRIPT Yancopoulos, G.D., Glass, D.J., 2004. The IGF-1/PI3K/Akt pathway prevents expression of muscle atrophy-induced ubiquitin ligases by inhibiting FOXO transcription factors. Mol. Cell 14, 395-403. Sun, Y., Guo, C. Y., Wang, D. D., Li, X. F., Zhang, Y. H., You, X. X., Shi, Q., Hu, G. J., Fang, C., Lin, H.R., Zhang, Y., 2016a. Transcriptome analysis reveals the molecular mechanisms underlying growth superiority in a novel grouper hybrid (Epinephelus fuscogutatus♀ × E. lanceolatus♂). BMC Genetics 17, 24. Sun, Y., Huang, Y., Hu, G.J., Zhang, X.H., Ruan, Z.Q., Zhao, X.M., Guo, C.Y., Tang, Z.J., Li, X.F.,
PT
You, X.X., Lin, H.R., Zhang, Y., Shi, Q., 2016b. Comparative transcriptomic study of muscle provides new insights into the growth superiority of a novel grouper hybrid. Plos One 11, e0168802.
RI
Verma, E., Chakraborty, S., Tiwari, B., Mishra, A.K., 2018. Transcriptional regulation of acetyl CoA and lipid synthesis by PII protein in Synechococcus PCC 7942. J. Basic Microb. 58,
SC
187-197.
Wang, C.S., Hartsuck, J.A., 1993. Bile salt-activated lipase. A multiple function lipolyt ic enzyme. Biochim. Biophys. Acta. 1166, 1-19.
NU
Wang, S., Tang, C.C, Tao, M., Qin, Q.B., Zhang, C., Luo, K.K., Zhao, R.R., Wang, J., Ren, L., Xiao, J., Hu, F.Z., Zhou, R., Duan, W., Liu, S.J., 2019. Establishment and application of
MA
distant hybridization technology in fish. Sci. China Life Sci. 62, 22-45. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57-63.
Xiao, J., Kang, X.W., Xie, L.H., Qin, Q.B., He, Z.L., Hu, F.Z., Zhang, C., Zhao, R.R., Wang, J.,
ED
Luo, K.K., Liu, Y., Liu, S.J., 2014. The fertility of the hybrid lineage derived from female Megalobrama amblycephala × male Culter alburnus. Anim. Reprod. Sci. 151, 61-70.
EP T
Xie, C., Mao, X.Z., Huang, J.J., Ding, Y., Wu, J.M., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L.P., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316-W322. Xu, W.N., Qian, Y., Li, X.F., Li, J.Y., Li, P.F., Cai, D.S., Liu, W.B., 2017. Effects of dietary biotin on growth performance and fatty acids metabolism in blunt snout bream, Megalobrama
AC C
amblycephala fed with different lipid levels diets. Aquaculture 479, 790-797. Zhai, R.R., Feng, Y., Wang, H.M., Zhan, X.D., Shen, X.H., Wu, W.M., Zhang, Y.X., Chen, D.B., Dai, G.X., Yang, Z.L., Cao, L.Y., Cheng, S.H., 2013. Transcriptome analys is of rice root heterosis by RNA-Seq. Bmc Genomics 14, 19. Zhang, C., Lin, C., Fu, F., Zhong, X., Peng, B., Yan, H., Zhang, J., Zhang, W., Wang, P., Ding, X., Zhang, W., Ding, X., Zhang, W., Zhao, L., 2017. Comparative transcriptome analysis of flower heterosis in two soybean F1 hybrids by RNA-seq. PloS one 12, e0181061. Zhang, H. Y., He, H., Chen, L.B., Li, L., Liang, M.Z., Wang, X.F., Liu, X.G., He, G.M., Chen, R.S., Ma, L.G., Deng, X.W., 2008. A genome-wide transcription analysis reveals a close correlation of promoter INDEL polymorphism and heterotic gene expression in rice hybrids. Mol. Plant 1, 720-731. Zheng, G.D., Guo, D.D., Wu, C.B., Chen, J., Jiang, X.Y., Zou, S.M., 2019. The obvious heterosis and genetic characters of intergeneric cross and backcross juveniles between blunt snout
23
ACCEPTED MANUSCRIPT bream (Megalobrama amblycephala) and topmouth culter (Culter alburnus). Aquac. Res. 50, 1634-1643. Zheng, G.D., Wang, C.L., Guo, D.D., Jiang, X.Y., Zou, S.M., 2017. Ploidy level and performance in meiotic gynogenetic offsprings of grass carp using UV-irradiated blunt snout bream sperm. Aquaculture and Fisheries 2, 213-219. Zheng, G.D., Zhang, Q.Q., Li, F.G., Chen, J., Jiang, X.Y., Zou, S.M., 2015. Genetic characteris tics and growth performance of different Megalobrama amblycephala (♀) × Erythroculter
NU
SC
RI
PT
ilishaeformis (♂) hybrids. J. Fish. Sci. China 22, 402-409.
MA
Figure legends Fig. 1. Landscape of unigene distribution in the backcross BC. (A) E-value distribution of unigenes searched against public databases with an E-value cut-off of 1E-5. (B) Identity distribution of unigenes searched against public databases with an E-value cut-off of 1E-5. (C) Unigenes conserved in BC and model species. Unigenes of the BC were characterized by species by using BLASTX.
ED
searching against public databases. The number of BC homologous genes identified in other
EP T
Fig. 2. DEGs analysis and volcano plot for ‘BC vs. MA’ (A) and ‘BC vs. CA’ (B). The x-axis is the value of Log2 (Fold Change), and the y-axis is the value of Log2 (p-value). The red and yellow dots reveal the up-regulated DEGs, the blue and light blue dots reveal the down-regulated DEGs. MA, CA and BC denote M. amblycephala, C. alburnus and their backcrossing offsprings,
AC C
respectively.
Fig. 3. DEGs and the overlaps of DEGs between ‘BC vs. MA’ and ‘BC vs. CA’. Fig. 4. GO enrichment of DEGs in hepatopancreas of BC and their parents. (A) GO enrichment of DEGs of ‘BC vs. MA’. (B) GO enrichment of DEGs of ‘BC vs. CA’. Asterisk showed the most enriched GO terms in each category. Fig. 5. Hierarchical cluster analysis of DEGs involved in the growth superiority. The color key represents FPKM normalized log2 transformed counts in hepatopancreas of three fish species. Changes in expression levels are shown using color scales with saturation at > 2-fold changes. Red and green gradients indicate a increase and decrease of DEGs, respectively.
24
ACCEPTED MANUSCRIPT Fig. 6. A summary of the SSRs identified from the MA, CA and BC. Fig. 7. Distribution of putative single nucleotide polymorphisms (SNP) in the MA, CA and BC. Fig. 8. Validation of RNA-seq data of ten DEGs by qRT-PCR. β-actin was used as an internal control and each value represents average of three separate biological replicates. Fig. 9. The overview map of molecular mechanisms underlying growth superiority of the BC.
PT
Blue arrows denote the DEGs of ‘BC vs. MA’. Red arrows denote the DEGs ‘BC vs. CA’. Up or down arrows stand for up- or down-regulation in the BC compared with their parent. CM: cell membrane; MM: mitochondrial membrane; ERM: endoplasmic reticulum membrane; ACC: acetyl
RI
CoA carboxylase; oaa: oxaloacetate.
SC
Supplementary Information Supplementary Fig. 1. Length distribution of transcripts and unigenes in BC (A, D), MA (B, E)
NU
and CA (C, F) respectively.
Supplementary Fig. 2. COG functional classification of unigenes in BC. A total of 7,406
MA
assembled unigenes were annotated and assigned to 25 functional categories. Supplementary Fig. 3. KEGG pathway enrichment analysis of the DEGs of ‘BC vs. MA’. The x-axis is the significant enriched KEGG pathway classification and y-axis is P-value of
ED
enrichment.
enrichment.
EP T
Supplementary Fig. 4. KEGG pathway enrichment analysis of the DEGs of ‘BC vs. CA’. The x-axis is the significant enriched KEGG pathway classification and y-axis is P-value of
Supplementary Table 1. COG, KOG and NOG functional classification of unigenes in BC, MA
AC C
and CA, respectively.
25
ACCEPTED MANUSCRIPT Table 1. Genes and specific primers used for quantitative real-time PCR. Forward primer sequence (5’-3’)
Reverse primer sequence (5’-3’)
TRY
GGGTGTCTGATGACCTTAGTG
TGTCTGCTGCTCATTGCT
ELA1
TCTCCTCCTCCTTCAGGTTATG
CTTGTGGAGGCAGCCTTATT
CPA2
GCTTTCACTCACACCAATAACC
CCAGCATCCCAGTTCCTATT
BAL
GAGGAGATCGCAAAGAAGGTAG
ATCCACATTCCCAGCAAGAG
IGF1
CAAGAGAGGAGTTTGCAGTGA
GCAAGGGTTCCATCTGGTATAA
IGFBP2b
GCTGACAGAGGGCTAGATTTG
CTGCAGGCAAACTTGTGTTTAT
RAPTOR
CTGCCTCACTACACCAATCAA
GCCTGGGATCTTCTCAATCAA
EIF4E
CAGATGGGCTCTCTGGTATTTC
AGTTTACTGGGCTGCTGTATG
FAS
AGGTGTCCGAGAATGGAAATC
ELOVL1
CCGGCGTACAGTACAAAGAA
AGACCCAAGGTTGAAGGATTAC
β-actin
CGTGCTGTTTTCCCTTCCATT
CAATACCGTGCTCAAAGGATACTT
SC
RI
PT
Gene name
AC C
EP T
ED
MA
NU
GTCTGCAGGCATTGGTTTATTC
26
ACCEPTED MANUSCRIPT
Table 2. Growth performance and digestive enzyme activities (U g-1 tissue protein) Group
IW
FW
SGR
WGR
HSI
CF
Protease
Lipase
BC
6.50±0.35
50.70±3.47a
2.28a
680.00a
1.48±0.17ab
2.10±0.02a
31.60±1.34a
22.18±2.02a
MA
6.50±0.34
44.81±3.24b
2.14ab
589.38b
1.68±3.24a
1.86±0.01ab
23.50±1.10b
20.16±1.50b
CA
6.50±0.62
28.25±2.02c
1.63c
334.61c
1.07±0.11b
1.11±0.03c
25.43±1.15b
18.60±2.10b
PT
Notes: IW (g): initial body weight, FW (g): final body weight, Weight gain rate (WGR, %) = (FW − IW)×100 / IW, Specific growth rate (SGR, %) = (lnFW − lnIW)
RI
×100/number of culture days, Hepatosomatic index (HSI, %) = 100×hepatopancreas
SC
weight/body weight, Condition factor (CF, %) = 100×(body weight, g)/(body length, cm)3 . Digestive enzyme activities were expressed as micromoles of hydrolyzed
NU
substrate min-1 g-1 tissue protein (U g-1 tissue protein). Means in the same column with
AC C
EP T
ED
MA
different superscripts are significantly different (P < 0.05).
27
ACCEPTED MANUSCRIPT Table 3. Characterization of raw/clean data and mapping rate Total raw
# of clean
Total clean
Number of reads
Mapping
reads
read size (bp)
reads
read size (bp)
mapped to assembly
ratea (%)
BC
50,064,120
7,559,682,120
49,063,830
7,316,679,170
40,859,480
83.28
MA
56,209,354
8,487,612,454
55,225,738
8,244,390,148
49,405,514
89.46
CA
52,814,668
7,975,014,868
51,725,452
7,714,644,875
Average
53,029,381
8,007,436,481
52,005,007
7,758,571,398
Total
159,088,142
24,022,309,442
156,015,020 23,275,714,193
RI
SC
46,230,046
89.38
45,498,347
87.37
136,495,040
——
Reads with ≤ 5 base mismatches were counted when mapped to the reference
NU
a
PT
# of raw
EP T
ED
MA
sequences.
AC C
Sample
28
ACCEPTED MANUSCRIPT Table 4. Summary of cDNA sequences of the BC and their parents
MA
CA
BC
MA
CA
Total sequence No.
77102
63367
67087
62248
55182
57366
Total length (bp)
71136472
52919414
56482685
50596979
41110483
42828640
Largest length (bp)
16621
15752
15127
16621
15752
15127
Smallest length (bp)
201
201
201
201
201
201
Average length (bp)
923
835
842
813
745
747
N50
1799
1524
1601
1581
1288
1349
GC%
45.43
45.69
45.57
45.34
45.48
45.36
AC C
EP T
ED
MA
NU
SC
BC
PT
Unigenes
RI
Transcripts
29
ACCEPTED MANUSCRIPT
Gene name
BC-FPKM
MA-FPKM
CA-FKPM
Digestive
TRY
4614.88±65.81
411.07±6.83
1658.76±35.34
enzyme
ELA1
1327.78±25.27
145.72±5.47
400.73±11.26
CTRL1
564.14±18.29
6.51±0.84
125.65±11.31
CPA2
92.67±5.87
5.03±0.18
20.96±2.45
BAL
425.73±20.28
20.12±2.37
61.82±8.46
IGF1
28.02±1.05
2.73±0.84
3.11±0.93
IGF2a
82.29±5.86
27.62±3.63
13.79±2.86
IGFBP1
28.59±1.04
137.95±3.86
241.57±10.81
IGFBP2b
57.74±3.17
11.62±0.89
9.01±1.04
Protein
PI3KR
56.03±5.61
5.46±1.15
6.23±0.86
synthesis
RAPTOR
182.19±25.88
61.33±6.43
90.33±8.89
EIF4E
152.65±12.31
62.36±9.39
13.88±2.84
Fatty acid
CS
101.67±16.81
19.13±1.87
27.50±1.03
synthesis
MDH
979.76±32.12
404.63±15.23 469.19±21.13
1379.37±59.11
658.22±37.27 82.96±6.82
ELOVL1
33.86±5.45
14.14±4.72
6.65±1.37
ELOVL5
283.73±27.80
73.77±9.17
83.06±8.82
1213.94±55.88
604.20±30.28 35.71±3.91
AC C
EP T
FAS
ELOVL6
RI
SC
NU
ED
GH/IGF axis
PT
Function
MA
Table 5. FPKMs of digestive enzyme, IGF system and protein/fatty acid synthesis related genes in BC, MA and CA.
30
ACCEPTED MANUSCRIPT Highlights Intergeneric backcross BC with growth superiority and enhanced digestive enzyme activity were obtained between Megalobrama amblycephala and Culter alburnus.
PT
134 differentially expressed genes (DEGs) were identified between BC and
RI
their parents in hepatopancreas by transcriptome analysis.
SC
The DEGs related to the IGF system, digestive enzyme, protein/fatty acid
AC C
EP T
ED
MA
NU
synthesis may contribute to growth advantage of backcross BC.
31