Transcriptome Profiling and Molecular Pathway Analysis of Genes in Association with Salinity Adaptation in Nile Tilapia Oreochromis niloticus

Nile tilapia Oreochromis niloticus is a freshwater fish but can tolerate a wide range of salinities. The mechanism of salinity adaptation at the molecular level was studied using RNA-Seq to explore the molecular pathways in fish exposed to 0, 8, or 16 (practical salinity unit, psu). Based on the change of gene expressions, the differential genes unions from freshwater to saline water were classified into three categories. In the constant change category (1), steroid biosynthesis, steroid hormone biosynthesis, fat digestion and absorption, complement and coagulation cascades were significantly affected by salinity indicating the pivotal roles of sterol-related pathways in response to salinity stress. In the change-then-stable category (2), ribosomes, oxidative phosphorylation, signaling pathways for peroxisome proliferator activated receptors, and fat digestion and absorption changed significantly with increasing salinity, showing sensitivity to salinity variation in the environment and a responding threshold to salinity change. In the stable-then-change category (3), protein export, protein processing in endoplasmic reticulum, tight junction, thyroid hormone synthesis, antigen processing and presentation, glycolysis/gluconeogenesis and glycosaminoglycan biosynthesis—keratan sulfate were the significantly changed pathways, suggesting that these pathways were less sensitive to salinity variation. This study reveals fundamental mechanism of the molecular response to salinity adaptation in O. niloticus, and provides a general guidance to understand saline acclimation in O. niloticus.


Introduction
Salinity is one of the most significant factors regulating distribution, abundance and diversity of aquatic animals [1]. Due to diverse distributions of aquatic animals from freshwater to brackish or marine water, various physiological strategies have been evoluted for salinity adaptation [2,3]. To reveal the adaptive processes of aquatic animals in salinity acclimatization, most studies have focused on growth [3], distribution [4], osmoregulation [5], production [6] and physiological responses [7,8]. However, the underlying mechanism of salinity adaption in fish has not been well understood especially at the molecular level such as integrated molecular pathway responses. The existing literature in molecular biology is limited to specific gene cloning, function determination and biological pathways [9,10], and little is known on the overall responsive pathways relevant to adaptive mechanism in fish to salinity changes.
Oreochromis niloticus is a unique model species to study salinity adaptation as it can live in a wide range of salinities [11,12]. In the past decade, research on O. niloticus in brackish water has been confined to the influence of salinity on physiology, development and breeding [13][14][15]. Recently, studies in molecular biology have examined the influence of salinity on the expression of target genes such as the mRNA expression of Na + , K + -ATPase [16], renin-angiotensin system genes [17], transient receptor potential vanilloid 4 [18], growth hormone and somatolactin [19]. These fragmental studies have provided a basis to further explore adaptable strategies of O. niloticus to saline water. Therefore, there is a need to further study saline acclimatization of O. niloticus at a transcriptional level to reveal more fundamental mechanism in osmoregulation.
In fish, hepatopancreas is an important organ in energy metabolism [20] and detoxification [21,22]. A recent study shows that hepatopancreas can help maintain salt and fluid balance under salinity chanlenge in Senegalese sole (Solea senegalensis) [23]. Salinity challenge adds more stress on aquatic animals [3], and therefore more energy is needed to maintain homeostasis [24]. As hepatopancreas plays mutiple functions in fish, it is an ideal organ to test its reponse to salinty challenge through a comprehensive molecular approach.
With the rapid advances in molecular technologies, it has become possible to explore the ecological and physiological mechanisms regulating distribution and function of aquatic organisms [25]. Multiple approaches using transcriptome, digital gene expression and proteomics have been developed to further understand the molecular change in cells and tissues of fish under salinity challenge [26,27]. The transcription profiling method has been used to investigate the change of reference genes under environmental stress [26,28]. The transcriptomic analysis is essential to explain the underlying functional elements of the genome and can lay a foundation to help understand the ability of an organism to cope with various stress [29]. With the emergence of transcriptome sequencing, RNA-Seq has significantly improved the gene coverage and increased the sensitivity for differentially expressed genes [30].
Therefore, the present study used RNA-Seq to reveal the hepatopancreas transcriptome differences of O. niloticus at different salinities. The pathways and genes that respond to salinity were obtained and analyzed. The results will provide an insight into the underlying mechanism of salinity acclimation in O. niloticus and its homologous species.  (Table 1).

Differentially expressed genes and series clusters
The identification of genes was based on the Nile tilapia genome. To provide the gene annotation, the nucleic acid sequences of the genes were compared to the genomes of zebrafish, human, mouse and rat based on the amiGo database. Finally, a total of 296,000 genes were annotated. A total of 934, 1087 and 734 genes were differently expressed with either a fold change >2 or a fold change <0.5 (P < 0.05, FDR < 0.05) in the 0 vs 8 psu set, 0 vs16 psu set and 8 vs 16 psu set, respectively. Based on the comparison among these three groups, we obtained 1852 genes in differential gene unions and conducted eight types of unique model expression tendencies according to the amount of mRNA in relevant genes (Fig 1A and 1B). Tendencies 1, 6 and 7 were significantly different as calculated with Fisher's exact test and the multiple comparison tests (P < 0.05). Tendency 1 contained the most differentially expressed genes. Because tendency 2 and tendency 5 shared the same expressed genes without difference between the freshwater control and 16 psu, the analysis of these two tendencies was meaningless and not conducted.

GO and Pathway analysis
Gene ontology (GO) analysis involved annotated genes from functionally known species, and the gene products were divided into three categories: molecular function, cellular component and biological process. The unions of significantly enriched GO terms (P < 0.05) under biological processes containing three groups were divided according to different tendencies. There were 229 GO terms in tendency 0, 210 terms in tendency 1, 145 terms in tendency 3, 265 terms in tendency 4, 122 terms in tendency 6, and 152 terms in tendency 7 ( Table 2).
The Kyoto encyclopedia of genes and genomes database was used to obtain significantly changed pathways containing differentially expressed genes. We divided the 6 tendencies into 3 categories (Tables 3-5): constant change (containing tendencies 0 and 7), change-then-stable (containing tendencies 1 and 6) and stable-then-change (containing tendencies 3 and 4). Because tendencies 2 and 5 had no difference between the freshwater control and 16 psu, these pathways were not further analyzed. In the constant-change category, the complement and coagulation cascades contained four genes; the steroid hormone biosynthesis, steroid biosynthesis, and ovarian steroidogenesis contained 12 genes; and fat digestion and absorption, vitamin digestion and absorption and retinol metabolism contained eight genes. The stable-thenchange category contained 31 pathways such as the biosynthesis of unsaturated fatty acids, fatty acid elongation, protein processing in endoplasmic reticulum, glycolysis/gluconeogenesis, pyruvate metabolism and tight junction. The change-then-stable category involved 29 pathways related to lipid metabolism, cell cycle and oxidative phosphorylation.

Gene-act network and co-expression network analysis
A total of 162 differentially expressed genes were involved in the constant change category, and they were used to build a gene-act network profile. A total of 72 genes were down-regulated, while 90 genes were up-regulated. The main subnetworks included fat digestion and absorption, glycolysis/gluconeogenesis, steroid biosynthesis, complement and coagulation cascades, endoplasmic reticulum activity, cell connection and signal transport (Fig 2A-2G). The gene-act network profile of the change-then-stable category was composed of 113 differentially expressed genes with 36 up-regulated genes and 77 down-regulated genes. The major subnetworks contained lipid and glycerophospholipid metabolism, glucose utilization, protein and amino acid metabolism, PI3K signaling pathways, nucleotide metabolism, oxidative stress and cytoskeleton (Fig 3A-3G). A total of 101 differentially expressed genes composed the stable-then-change category, with 20 down-regulated genes and 81 up-regulated genes. The chief sub-networks which contained in this category were the signaling pathway of PI3K and p53, steroid hormone synthesis and oxidative stress, fat synthesis and glycerophospholipid metabolism, cytoskeleton, endoplasmic reticulum activity, sugar utilization and pyruvate metabolism (Fig 4A-4F). The gene co-expression profile of the constant-change category contained 293 genes that were differentially expressed relative to others. To locate the core regulatory genes that were involved in the molecular response to salinity, core effect factors were determined by the degree differences between the freshwater and 16 psu (Fig 5). There were 40 genes in the red ball with 8 in the k-core, indicating that these genes would have the broadest contact with the surrounding genes and share the same expression tendencies under the same circumstances.

Experimental validation
Twenty five differentially expressed genes were randomly selected to validate the reliability of the RNA transcriptome in this study. Melting-curve analysis revealed a single product for all of the tested genes. Log 2 FCs from qPCR were compared with the results of RNA-Seq expression analysis (Fig 6). The results of Q-PCR and RNA-Seq showed a correlation coefficient of 0.93, indicating the credible RNA-Seq results.

Pathway analysis in the constant-change category
The pathways in the constant-change category are sensitive to salinity fluctuation. When the pathways respond positively to salinity fluctuation, these intermediate products are required for the basic physiological activities [31]. Therefore, the stability maintenance of these pathways is primarily important in salinity acclimation of tilapia. Due to the sensitivity of these pathways, they may play critical roles in response to a variable environment, and it is    [32,33]. The steroid regulatory metabolism under salinity stress in O. niloticus is discussed below.
Biosynthesis of cholesterol is centered in the relationship between osmoregulation and steroid as the cholesterol induced by salinity is related to the physical properties of cell membranes [34] and steroid hormone in ovarian development. Steroid hormone demonstrates the osmoregulation ability in red blood cells against hypotonic hemolysis in dog [35,36]. In the present study, in a saline environment, the ovarian steroidogenesis pathway activated the cAMP signal pathway and stimulated adenylate cyclase for cAMP-regulated gene transcription (S2 Fig) in tilapia hepatopancreas [37]. The cAMP signaling pathway stimulates the production of arachidonic acid metabolites to regulate hormone production, including cortisol and glucagon, osmoregulation and cellular fatty acid signaling in tilapia [37,38].
The expression of gene sulfotransferase 2B1 promotes the production of cholesterol sulfate and dehydroepiandrosterone (DHEA) [39]. Previous research suggests that the organisms can utilize cholesterol sulfate to support platelet adhesion [40], regulate serine proteases to alter epidermal cell adhesion [41] and transmit with signal transduction by modulating selective protein kinase C isoforms and phosphatidylinositol 3-kinase (PI3K) [42,43]. The changes in the PI3K signal pathway and cell adhesion in this study suggest that cholesterol sulfate affects salinity domestication in tilapia. Research on the relationship between DHEA and salinity adaptation is scarce. The DHEA could be indirectly affected by salinity domestication based on the present research. In other animal models, DHEA and its metabolites constitute α subunit activators, which is one class of peroxisome proliferators-activated receptors (PPAR) found in nature [44,45] and also in the present study. The functions of DHEA, such as anti-inflammatory properties during inflammatory responses in mice [46], may be involved in the regulation of immune-related pathways. These findings could provide new research directions for studying DHEA in aquatic animals.
Combining the literature evidence with the RNA-Seq results in the present study, the sulfur metabolism response to salt adaptation is surmised nonspecifically by sulfonating pregnenolone, cholesterol and DHEA for efficiently regulating metabolism and steroids [47,48].
Lipid metabolism-related pathways. Under salinity stress, fish must have appropriate mechanisms to regulate osmotic balance and to maintain basic physiological functions [49,50] which has been verified in this study (S1 Table). It is reported that euryhaline fish and crab Chasmagnathus granulata can utilize lipids as an energy source when encountering osmotic stress [51,52], which is similar to the results of the present study. In this experiment, O. niloticus might decrease triglyceride production and reserve lysophosphatidic acid, monoacylglycerols and fatty acids to spare energy for osmoregulation which is alike to the results of other study [53].
Immune-related pathways. Retinol metabolism, vitamin digestion and absorption, porphyrin and chlorophyll metabolism and the metabolism of xenobiotics by cytochrome P450 were gradually down-regulated while the complement and coagulation cascades were up-regulated by salt pressure in this study. The immune-related pathways in this category mainly cover the requirements of antioxidant and complement reactions that are altered by appropriate timing to salinity variation. The carotenoids (including β-carotene), retinol (vitamin A) and αtocopherol (vitamin E), which are lipid-soluble antioxidants in aquatic animals [54], are utilized to contribute to antioxidants under salinity stress in the hepatopancreas. The transcripts of two genes encoding proteins involved in porphyrin and chlorophyll metabolism are changed and they can participate in oxidation-reduction reactions in freshwater fish [55]. In fish, the cytochrome P450 monooxygenase system metabolizes a large number of xenobiotic compounds to perform detoxification [56]. In addition, plasminogen, complement 3 and complement 6 [57] are highly expressed with salinity elevation to cope with anti-physiological challenges in the O. niloticus hepatopancreas.
Osmoregulation-related pathways. Constant salinity change provokes transcriptional up-regulation of genes involved in arginine and proline metabolism and protein digestion and absorption. In aquatic insect larvae, free amino acids such as arginine and proline act as regulators for osmoregulation [58]. In blue crab Callinectes sapidus, amino acids participate in cell volume regulation. However, the ways how these amino acids work in osmoregulation are not clear [59]. The arginine and proline metabolisms are also involved in transcriptome analysis of Chinese mitten crab under salinity stress [26]. Therefore, it is reasonable to deduce that constant salinity change provokes arginine and proline metabolisms which are helpful in conservative regulation under salinity stress in tilapia and other species.

Pathway analysis in the change-then-stable category
The change-then-stable category displays pathways that are changed in tilapia transcriptome and that are sensitive to salt stress, but the changes would stop when the salinity concentration exceeds a threshold value. The reasons for maintaining stable stages in these pathways are possibly that (1) the compensation approaches can cover the environmental changes; and (2) the management of these pathways has reached a limit, implying that the condition may be stressful to the organism.
Cell activity-related pathway. The ribosome is most significantly changed in tendency 1 in this study. Ribosomal proteins are ubiquitous and abundant to bind to RNA, and participate in balancing the synthesis of the RNA and ribosome protein components [60]. The relationships between osmoregulation and ribosome in animals are little known, but the ribosome in plants can decrease under osmotic stress [61,62] which is consistent with the present study and provides evolutionary evidence on the similarity between plants and animals for the ribosome in regulating osmotic stress. Changes in free ribosomal proteins in ribosome composition could have a relationship with the p53 system to regulate physiological activities [63], which also agrees with the result of the present study in the stable-then-change category. The cellular responses to osmotic stress affect the cell cycle and closely link to the DNA conformation and DNA activity [64]. Research on osmoregulation in chondrocytes has proven the reason for the change in the shape of the chondrocytes in the collagen fibrils [65]. This finding is consistent with the up-regulation of the collagen gene in the present study and shows that collagen participates in the salinity stress responses of the cells and tissues. These pathways could be involved in the regulatory mechanisms when an organism is under salt stress.
Immune-related pathways. Immune-related pathways are related to proteasome and oxidative phosphorylation with the characteristics of the change-then-stable category. In cells, the proteasome is inhibited by the p38 MAPK-dependent phosphorylation response to osmotic stress [66], which explains the down-regulation of proteasomes to alter the intracellular antigen process in the cytolytic immune response [67] under salt stress. The extracellular damage and change in osmotic pressure lead to an increase in the production of intracellular reactive oxygen species [68]. Therefore, extra reactive oxygen species can switch on as an immune response against homeostasis disorder.  [69] to enable the reduction of fatty acid biosynthesis due to a lack of materials in the substrate. In response to fatty acid loss, steroid biosynthesis and the PPAR signaling pathway may be down-regulated to decrease apolipoprotein and the expression of long-chain acyl-CoA synthetase to activate fatty acid β oxidation [70]. The indirect evidence of the PPAR signaling pathway mediated by osmotic stress is that angiotensin II is a hormone response to osmotic stress as well as its association with the down-regulation of PPARs [71]. Ether lipid is an emerging class of lipids and its function in osmoregulation is an enigma, though they act a platelet-activating factor in the membrane components in the brain and testis [72].
Pathway analysis in the stable-then-change category Transcriptome differences among 8 psu, 16 psu and the freshwater environment are not sufficiently stressful to activate the pathways that are involved in the stable-then-change category. These regulatory pathways can adapt to a low salinity environment when the salinity reaches a threshold that these paths cannot withstand the change and result in functional adaptation at the transcriptional level [73,74]. niloticus needs to upregulate some fatty acids synthesis pathways to meet the demands for these fatty acids. Similarly, the high level of dietary n-3 highly unsaturated fatty acids improves the tolerance of Chinese mitten crab larvae to salinity stress [75]. It is well known that fish are liable to congenital diabetes due to a low ability of carbohydrate utilization [76]. During seawater acclimation in Oreochromis mossambicus, the glycogen content in the hepatopancreas is decreased significantly after transfer to seawater, suggesting that salt acclimation promotes carbohydrate utilization [77], which supports the findings in the current study. Thus, it is suggested that lipids are the main source of energy supply under salinity domestication in O. niloticus.
Protein metabolism-related pathways. Protein pathways, including thyroid hormone synthesis, protein processing in endoplasmic reticulum, arginine and proline metabolism, protein export and the biosynthesis of amino acids belong to tendency 4. In fish, thyroid hormone is an amino acid derivative that acts as a regulator for metabolism, osmoregulation and salinity adaptation [78,79]. The thyroid hormone receptor regulates cholesterol and carbohydrate metabolism through direct actions on gene expression [80] as well as cross-talk with other nuclear receptors, including PPAR and the liver X receptor [81]. The result of this study suggests that material synthesis and transmission to counteract stress require a large amount of energy and related organelle activities.
Cell activity-related pathways. Tight junction and antigen processing and presentation are pathways of tendency 4, while other types of O-glycan biosynthesis and glycosaminoglycan biosynthesis-keratan sulfate belong to tendency 3. Long-term salinity acclimation leads to the increase of gene transcription in basic components of cytoskeleton such as tilapia claudins, actin and myosin, and the cytoskeleton can further interact with intercellular tight junctions [82]. The production of glycan in this study may play a role in information transmission due to its specificity in the immune response [83].
Signaling pathways. Signaling pathways such as the p53, PI3K-Akt and HIF-1 are increased by high salinity. P53 is a tumor suppressor protein and a critical component in cell cycle checkpoint response [84]. However, recent research suggests that the p53 level does not provoke an apparent increase in rainbow trout compared with mammals [85]. In aquatic animals, it is unclear whether p53 stimulated by salt adaptation is able to induce cancer cell proliferation. In this study, p53 transcription is not altered but the downstream genes are upregulated, and there is no sufficient evidence to prove that salinity-induced carcinogenesis in tilapia. Insulin-like growth factor 1 receptor transcription is up-regulated to generate the PI3K-Akt signaling pathway and HIF-1 signaling pathway [86]. However, these conclusions in other animal models are not verified in aquatic animals.

Conclusions
According to the survival and growth parameters of O. niloticus (S1 and S2 Tables), it is feasible to rear O. niloticus at saline water. From the molecular perspective as shown in Fig 7, during salinity acclimation, O. niloticus produced significant changes in amino acid metabolism and synthesis, oxidation, protein synthesis and degradation, energy material utilization, and signal transduction. Glycolysis and fatty acids are involved in the regulation of acetyl coenzyme A synthesis and metabolism to participate in the TCA cycle and produce ATP for energy supply. Acetyl coenzyme A also participates in cholesterol synthesis by adjusting the needs of ovarian steroids and steroid hormone synthesis. Ovarian steroidogenesis activates the cAMP signal pathway to regulate adenylate cyclase, downstream gene expression and arachidonic acid metabolites. Among these actions, adenylate cyclase catalyzes ATP into cAMP to support signal transmission and the downstream genes of cAMP signal pathway cover various physiological processes. Arachidonic acid metabolites play extensive roles in maintaining homeostasis. Steroid hormone biosynthesis produces cholesterol-containing DHEA and cholesterol sulfate, which in turn participate in the PPAR pathway, immune-related pathways, cell connections and the PI3K signal pathway. The synthesis and metabolism of some amino acids reflect the reaction of tilapia to maintain osmotic stability. Protein synthesis and the metabolism in organisms are a prerequisite before responding to an environmental salinity challenge. In short, the steroid hormones, osmoregulation, lipid metabolism and cell-connected components are critical measures for salinity domestication in aquatic animals.

Experimental design and sampling
Oreochromis niloticus juveniles (initial weight 0.88 ± 0.04 g) were purchased from a farm in Hainan, China. After acclimation in freshwater for two weeks, fish were stocked in nine freshwater tanks (66 × 63 × 40 cm) with 25 fish each. One set of three tanks with freshwater were used as control and other two sets of three tanks were acclimated to the salinities of 8 and 16 psu, respectively, by increasing the salinity at 4 psu per day using crystal sea salt. The salinity was daily measured with a salimeter (AZ Instrument Corp. Ltd, Co. AZ8371). Fish were exposed to three final salinities 0, 8, and 16 psu in triplicate for 8 weeks. During the experiment, the environment was maintained at 12:12 dark/light, 28 ± 1°C, 7.99 ± 0.23 pH and >6 mg/L dissolved oxygen. One third of the water volume was daily renewed with corresponding salinities. Tilapias were fed with a commercial diet (35% crude protein and 3% crude lipid, S3 Table) twice daily at 08:00 h and 16:00 h to apparent satiation. Feces and uneaten feed were daily removed with a siphon tube. At the end of the experiment, O. niloticus in all tanks were fasted for 24 h before sampling. All experiments were conducted under the standard code of protocol for the care and use of laboratory animals in China. This research project was approved by the Animal Ethics Committee of East China Normal University. Eight fish were randomly selected and were anesthetized in 30 ppm MS-222. The hepatopancreas of each fish was put in a sterile plastic tube with Trizol (Invitrogen) for RNA extraction.

RNA extraction
The total RNA of the hepatopancreas was extracted by the Trizol method (Invitrogen). The integrity of every RNA sample was tested using 1% agarose gel electrophoresis to ensure that the RNA was integral with three distinct and bright stripes. In each treatment, the eight qualified RNA samples from each treatment group were mixed with the same amount of RNA to make sure the veracities and universalities. The quality of mixed RNA was further checked on a Bioanalyzer 2200 (Agilent Technologies, Santa Clara, CA, USA) with a RIN e value of >8, which is acceptable for cDNA library construction (S10 Fig). The mRNA is purified using the Dynabeads mRNA Purification Kit (Life tech, Cat. no. 1264684).

cDNA library construction and sequencing
The RNA-Seq was conducted by Novel Bioinformatics Co., Ltd with Sanger / Illumina 1.9. The cDNA library was prepared using the Ion Total RNA-Seq Kit v2 (Life technologies, Cat. no. 4479789) with 5 μg total RNA following the manufacturer's instructions. The mRNA with poly (A) was isolated with Dynabeads (Life technologies, USA), fragmented with RNaseIII and then purified. The fragmented RNA was added and ligated with ion adaptor. Double-stranded cDNA was synthesized and purified by the magnetic bead based method. The molar concentration of the purified cDNA was detected for each cDNA library. The filtering of the raw reads was produced with FAST-QC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Samples were processed on a OneTouch 2 instrument and enriched on a OneTouch 2 ES station for preparing the template-positive Ion PI

Quality control and mapping
The raw sequencing data were also evaluated by FAST-QC. The evaluation metrics scrutinized the nature of the data before subsequent variant evaluation. MapSplice software (MapSplice v2.1.8) was used, which is a tool for RNA-seq data analysis on mapping the RNA-seq data to the genome of O. niloticus Orenil1.1. MapSplice is the core program in Bowtie 1.0.0 that can quickly identify splicing and cutting the exon-exon [87]. The accession number of RNA-Seq is GSE67671.

Differential expression and series cluster
Different gene expression was screened by the EB-Seq algorithm [88]. The differential gene screening criteria was fold change >2 or fold change <0.5 (P <0.05, and false positive rate (FDR) <0.05). If multiple transcripts existed in a gene, the longest transcript was selected to calculate the sequencing depth and expression to ensure the result accuracy. The whole differential gene union of the three sets was obtained. Any two of the three sets (0 vs 8, 0 vs 16 and 8 vs 16 psu) were compared for the signal values of genes in the differential gene union. A series of expression tendency models of the differential genes in three groups were obtained. The minimum correlation coefficient of the state tendency analysis is 0.85. Significant models had a higher probability than expected by Fisher's exact test and the multiple comparison tests [89].

GO analysis and pathway analysis
Gene Ontology (GO) annotation of differential genes was based on the AmiGO database (http://amigo.geneontology.org/amigo). Using Fisher's exact test and χ 2 test [89], we calculated the result of the multiple hypothesis testing correction to obtaining the FDR, which was produced to correct the P-value. The FDR was defined as FDR = 1-N K /T, where N K refers to the number of Fisher's test P-values less than the χ 2 test P-values. Thus, the differentially expressed genes in significantly enriched GO terms were selected at P-value <0.01 and FDR <0.05. Pathway annotations of differential genes were downloaded from Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/). A Fisher exact test was used to find the significant enrichment pathways. The resulting P-values were adjusted using the BH FDR algorithm [89]. Pathway categories with a P-value <0.05 were reported. Enrichment provided a measure of the significance of the function because more enrichment represents more significant pathways in the experiment.

Gene-act network and co-expression network
We used the KEGG database to build the network of genes according to the relationship among the genes, proteins and compounds in the database [90][91][92][93]. Then, the gene co-expression networks were presented to find the interactions among the genes [94]. The gene coexpression networks were built according to the normalized signal intensity of specific expression genes. For each pair of genes, we calculated the Pearson correlation (P-value <0.05, Pearson >0.999 or <-0.999) and chose the significant correlation pairs to construct the network [95]. Within the network analysis, degree centrality is the simplest and most important measure of the centrality of a gene within a network to determine the relative importance. The degree centrality is defined as the number of links that one node has to other [96]. Moreover, to study the properties of the networks, k-cores in the graph theory were introduced as a method of simplifying the graph topology analysis. A k-core of a protein-protein interaction network usually contains cohesive groups of proteins [97,98]. The purpose of network structure analysis was to locate the core regulatory factors (genes) that were in one network. Core regulatory factors connected most of the adjacent genes and had the largest degree. While considering different networks, the core regulatory factors were determined by the degree differences between two classes of samples. They usually had the largest degree of differences.

QPCR validation
The reverse transcription of the three groups of extracted RNA was followed by PrimeScript RT Master Mix (TAKARA, Code No. RR036A). The reaction volume was 20 μL at the end and consisted of 4 μL 5 × PrimeScript RT Master Mix (Perfect Real Time), 1 μg total mRNA, and a supplement of nuclease-free water. The protocol of reverse transcription was 37°C for 15 min, and 85°C for 5 sec. Twenty five genes were randomly selected to design specific primers with primer 6.0 ( Table 6). The validation of RNA-Seq was conducted in a final volume of 20 μL that contained 3 μL of cDNA, 0.4 μM of each primer, 10 μL of ULtraSYBR Mixture and 6.2 μL of nucleasefree water. PCR amplifications were performed by the Bio-Rad CFX96 RealTime PCR system (Bio-Rad, US). The reaction program was the following: 95°C for 30 sec, 50 cycles (94°C for 15 sec, 58°C for 20 sec), and 72°C for 20 sec. The dissolution curve temperature is 60.0-95.0°C, which was increased by 0.5°C per 0.05 sec. The cDNA of each group with three parallels was normalized with the EF1A gene [99]. Then, the 2 -ΔΔCt method processing data was used to obtain the fold change [100].
Supporting Information