Adult pompano (T. ovatus) (body weight, 300 ± 20 g, body length: 22 ± 3 cm) were purchased from the Dongfeng Market (Zhanjiang, Guangdong, China), female fish were selected according to the observed gonads. The sbGnRH peptide used for this study was purchased from GL Biochem (Shanghai, China). sbGnRH peptides were dissolved in physiological saline. All the fish were anesthetized with 100 mg/L tricaine methane sulfonate (MS-222, Sigma,St. Louis, MO, USA) and dissected. Animal experiments were conducted in accordance with the guidelines approved by the Animal Research and Ethics Committees of Fisheries College of Guangdong Ocean University, Zhanjiang, China.

The dose and sampling time of sbGnRH peptides injection were selected according to our previous research results (10). Female adult fish were intraperitoneally injected with 10 ng/g of sbGnRH peptides according to the average body weight, while the control group was intraperitoneally injected with physiological saline. Each group contained three fishes. Three hours after injection, the selected fish were anesthetized and the liver was collected, frozen in liquid nitrogen and later stored at -80°C.

Liver samples were accurately weighed and homogenized (diluted 1:10) in Phosphate buffered saline (PBS) (Solarbio, Beijing, China) with a tissue crusher. The homogenization was 900 × g centrifuge for 10 min at 4°C and retained the supernatant. The protein content of the homogenate was measured using Folin-phenol reagent (16, 17). The contents of albumin (ALB), cholesterol (CHO), high-density lipoprotein cholesterol (HDLC), low-density lipoprotein cholesterol (LDLC), glucose (GLUC), creatine kinase (CK), iron (IRON), magnesium (MG), calcium (CA), phosphate (PHOS), alkaline phosphatase (ALP), aspartate aminotransferase (AST), alanine aminotransferase (ALT) and creatinine (CREA) were detected by Roche Cobas C311 automatic biochemical analyzer (Shanghai, China). The kits used were Cobas C311 products.

According to the product instructions, amylase (AMS), glucose-6-phosphate dehydrogenase (G-6-PD), malic enzyme (ME), pyruvate kinase (PK), acyl CoA oxidase (ACO), lipase (LPS), fatty acid synthase (FAS), carnitine acyl carnitine translocation (CACT), acetyl CoA carboxylase (ACC), lipoprotein lipase (LPL), superoxide dismutase (SOD), Phospholipid hydroperoxide glutathione peroxidase (GSH-PX), catalase (CAT) and malondialdehyde (MDA) were detected. See Table S1 for product numbers.

Six liver samples (3 replicates per group) were used to prepare a transcriptome (RNA-seq) sequencing libraries by the Biomarker Technologies. The RNA-seq process was performed as follows: 1 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext^®Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. The library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. After which, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bio analyzer 2100 system. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform and paired-end reads were generated. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity (18) with the min_kmer_cov set to 2 by default and all other parameters set to default.

UniGene sequence was compared with Non-Redundant protein (NR) (19), Swiss prot (20), Clusters of Orthologous Group (COG) (21), Eukaryotic Orthologous Groups (KOG) (22) and evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG 4.5) (23) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (24) database alignment by using diamond (25) software. Kobas (26) was used to obtain the KEGG ontology results of UniGene in KEGG. Interproscan (27) used the database integrated by interpro to analyze the Gene Ontology (GO) (28) ontology results of new genes. After predicting the amino acid sequence of UniGene, HMMER (29) software was used to compare results obtained with Protein Families (Pfam) (30) database to obtain the annotation information of UniGene. Bowtie (31) was used to compare the sequenced reads with the UniGene library. According to the comparison results, the expression level was estimated in combination with RNA-Seq by Expectation-Maximization (RSEM) (32). The expression abundance of the corresponding UniGene was expressed by fragments per kilobase of exon model per million reads mapped (FPKM) value. In the process of differential expression analysis, the recognized and effective Benjamin Hochberg method is used to correct the significance p-value obtained from the original hypothesis test. The corrected p-value, FDR (false discovery rate), is then used as the key index of differential expression gene screening to reduce the false positive caused by independent statistical hypothesis test on the expression value of a large number of genes (33). In the screening process, p-value < 0.01 and the difference multiple FC (fold change) ≥ 1.5 are used as the screening criteria.

The expression pattern of differentially expressed genes (DEGs) in RNA-seq analysis was verified by RT-qPCR (34). Total RNA was extracted from liver tissues of the control group and sbGnRH peptides treated group using Trizol according to manufacturer’s instructions (Invitrogen, CA, USA). Each group had 3 repetitions. According to the PrimeScript RT Master Mix Perfect Real -Time Kit (Takara, China) manufacturer’s instructions, each sample has the same amount of RNA for reverse transcription (35, 36). RT-qPCR was performed with the LightCycle 480 system (Roche, Basel, Switzerland) using SYBR Premix Ex Taq II (TaKaRa Bio Inc., Shiga, Japan). Primer pairs used in this study were shown in Table 1 . The reference gene β-actin was used as an internal reference to normalize the mRNA level. The relative gene expression was calculated by 2^−ΔΔCt method.

Data were expressed as the mean ± standard error (SE). All statistical tests were performed using Statistical Package for the Social Sciences (SPSS) 19.0 (SPSS, Chicago, IL, USA). Significant differences in the data among groups were tested by one-way analysis of variance (ANOVA), followed by Duncan’s post-hoc test. The probability level lower than 0.05 (P < 0.05) indicated significance difference and lower than 0.01 (P < 0.01) is extremely significance.

The indexes of protein metabolism showed that there was no significant difference in the content of liver ALB in pompano injected with different concentrations of sbGnRH peptides (P > 0.05). Compared with the control group, the content of ALB in the treatment group increased ( Figure 1A ). In lipid metabolism, the content of CHO in the treatment group decreased significantly (P < 0.05), the content of HDLC increased significantly (P < 0.05), and the content of LDLC increased, but there was no significant difference between the treatment group and the control group (P > 0.05) ( Figure 1B ). In terms of energy metabolism, the contents of GLUC and CK in the treatment group increased significantly (P < 0.05) ( Figure 1C ). There was no significant difference in the contents of IRON, MG, CA and PHOS in the inorganic components (P > 0.05), but the contents of IRON, MG and PHOS increased and CA decreased in the treatment group ( Figure 1D ). The content of ALP in the liver function index increased, which had no significant difference with the control group (P > 0.05); The contents of AST and ALT were significantly higher than those in the control group (P < 0.05) ( Figure 1E ). The content of CREA in the renal function index increased significantly compared with the control group (P < 0.05) ( Figure 1F ).

After injection of sbGnRH peptides at the concentration of 10 ng/gbw for 3 h, the comparison results of digestive enzyme activities showed that the activities of AMS and LPS in the liver decreased, and there was no significant difference with the control group (P > 0.05) ( Figure 2A ). In terms of enzymes related to carbohydrate metabolism, the activities of PK and ACO decreased, but there was a significant difference in ACO compared with the control group (P < 0.05) ( Figure 2B ). In terms of antioxidant defense, GSH-PX activity decreased, and there was no significant difference with the control group (P > 0.05); SOD (P < 0.05) and CAT (P < 0.01) activities decreased significantly ( Figure 2C ). In terms of lipid metabolism, LPL activity decreased (P > 0.05), G-6-PD and FAS activities decreased significantly compared with the control group (P < 0.05), while HL activity increased (P > 0.05), CACT (P < 0.05), ME and ACC (P < 0.01) activities increased significantly compared with the control group ( Figure 2D ). In addition, a significant increase in MDA content reflecting oil rancidity was observed (P < 0.05) ( Figure 2E ).

Based on Sequencing by Synthesis (SBS) technology, the liver tissue of pompano was sequenced 3 hours after in vivo injection of 10 ng/gbw (treatment group) sbGnRH peptides using Illumina hiseq high-throughput sequencing platform. Raw reads were obtained from C-1, C-2, C-3, T-1, T-2, and T-3. After quality control and removal of some low-quality data, clean reads of 21,283,647、19,427,359、21,873,990、21,732,174、23,660,062 and 21,592,338 were obtained respectively. The sequencing quality results show that the Q20 value (percentage of sequences with a sequencing error rate less than 1%) of each group of samples exceeds 97%, the Q30 value (percentage of sequences with a sequencing error rate less than 0.1%) exceeds 92%, and the GC content was about 46.92% - 47.80% ( Table 2 ).

UniGene sequences were compared with COG, GO, KEGG, KOG, Pfam, Swiss-Prot, eggNOG and NR databases to obtain the annotation information of UniGene ( Table 3 ). 25,715 Unigene with annotation information were obtained. Among the results of all database alignment genes, 4527 (17.60%) were annotated in the COG database, 22280 (86.64%) were annotated in the GO database, 20570 (79.99%) were annotated in the KEGG database, 15706 (61.08%) in KOG database, 17050 (66.30%) in Pfam database, 11992 (46.63%) in Swiss-Prot database, 21702 (84.39%) in eggNOG database and 24797 (96.43%) in NR database.

UniGene sequences were aligned into the NR database. The species with the largest number were Seriola dumerili (7274, 29.33%), followed by Seriola lalandi (4534, 18.28%), Lates calcifer (2550, 10.28%), Danio rerio (1401, 5.65%), Echeneis naurates (812, 3.27%), Scophthalmus maximus (473, 1.91%), Larimichthys crocea (419, 1.69%), Morone saxatilis (300, 1.21%), Epilephelus lanceolatus (257, 1.04%), Stegastes partitus (231, 0.93%). This indicates that the species with the closest genetic relationship compared with the liver transcriptome sequencing of pompano is Seriola dumerili ( Figure 3 ).

In this study, 99 genes were screened and identified as differentially expressed genes, including 77 up-regulated genes and 22 down-regulated genes ( Figure 4A ). Hierarchical cluster analysis was performed on the screened differentially expressed genes. The results showed that the biological repeat individuals of liver tissue after treatment with different concentrations of sbGnRH peptides could gather respectively, indicating that the repeatability of the sample was more reliable ( Figure 4B ). Based on the expression amount of genes in different samples, the identified differentially expressed genes were annotated. The number of differentially expressed genes annotated was 88, including 21 in COG database, 74 in GO database, 74 in KEGG database, 50 in KOG database, 72 in Pfam database and 56 in Swiss-Prot database, eggNOG database annotation to 76, NR database annotation to 86 ( Table 4 ).

According to GO enrichment analysis, DEGs are divided into three main functional categories: cellular component (CC), molecular function (MF) and biological process (BP), including 15, 12 and 18 subcategories respectively ( Figure 5 ). Among them, the DEGs in the CC category were significantly enriched in cell (GO: 0005623), cell part (GO: 0044464) and organelle (GO: 0043226) GO terms. In the MF category, the DEGs were significantly enriched in catalytic activity (GO: 0003824) and binding (GO: 0005488) GO terms. The majority of DEGs in the BP category were associated with cellular processes (GO: 0009987), single-organism process (GO: 0044699) and biological regulation (GO: 0065007) GO terms.

The KEGG pathway enrichment analysis results of liver DEGs show that 41 pathways are enriched, of which the first 20 pathways with the most reliable enrichment significance (the smallest q-value) are shown in Figure 6A . These pathways and the typical genes involved can be divided into cellular processes, environmental information processing, genetic information processing, diseases, metabolism and organismal systems ( Figure 6B ). Among them, the following KEGG pathways related to disease immunity were significantly enriched: Dilated cardiomyopathy, Renal cell carcinoma, Salmonella infection and Pathways in cancer; the metabolism related KEGG pathway is significantly enriched: Fatty acid metabolism, Pyruvate metabolism, Glycosphingolipid biosynthesis - lacto and neolacto series, Fatty acid elongation, Carbon metabolism, Glyoxylate and dicarboxylate metabolism, Glycerophospholipid metabolism, Inositol phosphate metabolism, Ether lipid metabolism, Folate biosynthesis, Amino sugar and nucleotide sugar metabolism, Glycosaminoglycan biosynthesis - keratan sulfate, Oxidative phosphorylation and Retinol metabolism ( Figure 6B ). These functional classifications provide a theoretical basis for revealing the effects of sbGnRH on liver diseases, immune and metabolic functions of pompano.

In this experiment, six genes (vtg, gpr1, egr1, tep1, bcl9l and rbm34) were randomly screened for qPCR verification, in order to verify the accuracy of transcriptome sequencing results. The results showed that the expression of vtg, gpr1 and egr1 genes was up-regulated and the expression of tep1, bcl9l and rbm34 genes was down-regulated after sbGnRH peptides treatment. The qPCR results were consistent with the sequencing results, indicating the accuracy and specificity of transcriptome analysis ( Figure 7 ).