Deep Transcriptomic Analysis of Black Rockfish (Sebastes schlegelii) Provides New Insights on Responses to Acute Temperature Stress

In the present study, we conducted an RNA-Seq analysis to characterize the genes and pathways involved in acute thermal and cold stress responses in the liver of black rockfish, a viviparous teleost that has the ability to cope with a wide range of temperature changes. A total of 584 annotated differentially expressed genes (DEGs) were identified in all three comparisons (HT vs NT, HT vs LT and LT vs NT). Based on an enrichment analysis, DEGs with a potential role in stress accommodation were classified into several categories, including protein folding, metabolism, immune response, signal transduction, molecule transport, membrane, and cell proliferation/apoptosis. Considering that thermal stress has a greater effect than cold stress in black rockfish, 24 shared DEGs in the intersection of the HT vs LT and HT vs NT groups were enriched in 2 oxidation-related gene ontology (GO) terms. Nine important heat-stress-reducing pathways were significantly identified and classified into 3 classes: immune and infectious diseases, organismal immune system and endocrine system. Eight DEGs (early growth response protein 1, bile salt export pump, abcb11, hsp70a, rtp3, 1,25-dihydroxyvitamin d(3) 24-hydroxylase, apoa4, transcription factor jun-b-like and an uncharacterized gene) were observed among all three comparisons, strongly implying their potentially important roles in temperature stress responses.

Ecosystems are currently exposed to global warming and climate change. One of the most direct impacts of climate change on the marine ecosystem affects fisheries. It has been reported that the temperature of the upper ocean (0 to 700 m depth) has increased, rising with an average rate of 0.05 °C per decade since 1971. The rate of temperature change is highest near the surface of the ocean (>0.1 °C per decade in the upper 75 m from 1971 to 2010) 1 . Fish are poikilothermic aquatic animals whose body temperatures adapt to environmental temperatures to a certain degree, changes in water temperatures may affect their growth, survival, reproduction, development and physiological performances 2,3 .
The molecular mechanisms underlying temperature stress conditions have long been of interest. Temperature stress causes expression changes in a series of stress-responsive genes, such as genes regulating protein folding repair 4,5 , energy metabolism 6,7 , the oxidation reduction process 7 , and the control of the cell cycle 8,9 . The identification of stress-responsive genes and pathways is the first step to reveal the fundamental mechanisms of the response to thermal stress and to predict the capacity of fish to adapt to climate changes. The next-generation sequencing technology (NGS)-based RNA-Seq platform is considered to be a revolutionary and efficient tool for investigating stress-responsive genes, as it can quantify over millions of unknown transcripts at once. RNA-Seq has been applied in studies of responses to temperature stress in several fish species, such as catfish 7 , Australian rainbowfish 10 , and snow trout 11 . However, almost all of these studies focused on oviparous fish species.
Ovoviviparity is a unique fish reproduction mode, in which fertilized eggs cannot delivered from the female ovary until the embryos are mature. Black rockfish (Sebastes schlegelii), belonging to Scorpaenidae, is an Annotation and function analysis of liver transcripts. All transcripts were subjected to annotation analysis by a comparison with the Nr, Nt, KEGG, KO, Swiss-Prot, PFAM, GO and KOG database. The results in Table 2 show the number of annotated transcripts in each database. A total number of 109,302(50.38%) transcripts were annotated by at least one database, and 66,596 (50.38%) annotated transcripts showed a significant BLAST hit against Nr database.
For the 66,596 transcripts that matched against the Nr database, the most abundant BLAST hits were from fish species (35.5%) such as Larimichthys crocea (19%), Stegastes partitus (8.9%) and Notothenia coriiceps (7.6%), followed by some other species, Homo sapiens (11.3%), Schistosoma japonicum (8.7%)and others (44.5%) (Supplement 1a). The functional classification of transcriptome data is the primary requirement for the application of functional genomic approaches in fishery research. GO and KEGG analyses are currently the most popular methods used for the functional classification of transcriptomic sequences. Our results showed that Blast2Go assigned 47,427 transcripts 56 functional GO terms (Supplement 1b). Regarding the three primary ontology categories, BP represents the majority (25 terms) of annotations, followed by CC (20 terms) and MF (11 terms). Based on the analysis of level 2 GO terms, the GO terms in BP with the highest numbers of annotations were cellular process (GO:0009987), metabolic process (GO:0008152), single-organism process (GO:0044699), biological regulation (GO:0065007) and regulation of biological process (GO:0050789). For CC, cell (GO:0005623), cell part (GO:0044464), organelle (GO:0043226) and macromolecular complex (GO:0032991) contained the highest numbers of annotations. The GO terms related to MF with the highest number of annotations were binding (GO:0005488), catalytic activity (GO:0003824) and transporter activity (GO:0005215). KEGG analysis was performed to understand the higher order functional information of biological system 17 . Based on the analysis, a total of 34,140 sequences were annotated with five categories on 232 KEGG pathways (Supplement 1c).

Analysis of differentially expressed genes (DEGs). Identification of DEGs.
A total of 584 annotated transcripts showed significant differential expression in at least one comparison (HT vs NT, HT vs LT and LT vs NT) (adjusted q-value < 0.05). Among them, 362 (223 up-regulated and 139 down-regulated) differentially expressed genes (DEGs) were identified in the HT vs NT group, and 421 (198 up-regulated and 223 down-regulated) DEGs and 113 (74 up-regulated and 39 down-regulated) DEGs were found in the HT vs LT groups and the LT vs NT groups, respectively (Supplement 2). The heat map presents the differently expressed transcripts and shows that LT and NT clustered in one group, indicating that cold stress causes fewer changes than heat stress in black rockfish (Supplement 3).
GO enrichment analysis 18 was performed on the 584 DEGs. In the HT vs LT group, a total of 151 DEGs has represented significantly enriched GO classifications, including 52 DEGs assigned oxidation-reduction process (BP, GO:0055114), 51 DEGs assigned oxidoreductase activity (MF, GO:0016491), 15 DEGs in oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (MF, GO:0016705), 11 DEGs assigned heme binding (MF, GO:0020037), 11 DEGs assigned iron ion binding (MF, GO:0005506) and 11 DEGs in tetrapyrrole binding (MF, GO:0046906) (Fig. 1a). In HT vs NT group, only 83 DEGs represented an enriched GO classification, including 41 DEGs in oxidation-reduction process (BP, GO:0055114) and 42 DEGs assigned oxidoreductase activity (MF, GO:0016491) (Fig. 1b). There was no DEG enrichment observed for the GO classification of the LT vs NT group.
DEGs were mapped to several specific pathways by the KEGG pathway analysis, which included 194, 189 and 75 KEGG pathways in the HT vs LT group, HT vs NT group and LT vs NT group, respectively. Here, we present 26 significantly enriched pathways (q-value < 0.05); of these, the HT vs LT group mainly contains 11 pathways, including influenza A (KO: 05164), legionellosis (KO: 05134) and estrogen signaling pathway (KO: 04915); the HT vs LT group mainly contains 14 pathways including antigen processing and presentation (KO:04612) and NOD-like receptor signaling pathway (KO: 04621), and the LT vs NT group contains the PPAR signaling pathway (KO: 03320) (Supplement 4).
Based on the KEGG pathway analysis 17,19,20 and manual literature searches, a total of 245 candidate genes associated with stress responses and adaptation were identified in the 584 annotated DEGs; these candidate genes were filtered and classified into 8 functional categories including protein folding, metabolism, immune response, signal transduction, molecule transport, membrane, and cell proliferation/apoptosis (Supplement 5). A total of 100 genes with a |log2(fold change)| > 2 were selected from the 245 candidate genes and are listed in Table 3. The imputed putative functions of these genes are covered in the Discussion.

Comparative expression analysis in HT vs LT group, HT vs NT group and LT vs NT group.
Based on the 245 DEGs mentioned above, 196, 163 and 37 annotated DEGs were obtained from the HT vs LT group, HT vs NT group and LT vs NT group, respectively (Supplement 5) (Fig. 2).
For the intersection analysis among the different groups comparisons, a total of 123 DEGs were shared between the HT vs LT and the HT vs NT groups, which presented the maximum numbers of shared DEGs among  (Table 4) were shared between the HT vs LT and the HT vs NT groups (Table 4). GO enrichment analysis suggested 26 genes were enriched in the oxidation-reduction process term (13 up and 13 down DEGs) in BP and 26 were enriched in the oxidoreductase activity term (13 up and 13 down DEGs) in MT, including 24 DEGs sharing both GO terms (Supplement 2). Nine KEGG pathways were selected from the intersection of the HT vs LT group and the HT vs NT group, which were enriched in 11 and 14 KEGG pathways, respectively (q-value < 0.05) (Supplement 4). Among these 9 KEGG pathways, four pathways, influenza A (Fig. 3), legionellosis, inflammatory bowel disease (IBD) and measles were related to immune and infectious diseases. Four pathways, NOD-like receptor signaling pathway (Fig. 4), osteoclast differentiation, plant-pathogen interaction IBD and antigen processing and presentation, were related to organismal systems especially the immune system, and the estrogen signaling pathway (Fig. 5) was related to the endocrine system. There were only 37 DEGs (25 up-regulated and 12 down-regulated) identified in the LT vs NT group, which further indicated cold stress may cause less changes than thermal stress.

Validation of RNA-Seq results by qRT-PCR.
To validate the RNA-Seq results, 11 DEGs were randomly selected for qRT-PCR analysis (Supplement 6). The results showed that the qRT-PCR expression trends of the selected genes were significantly correlated with the RNA-Seq results (R 2 : 0.882-0.911). Generally, the RNA-Seq results were confirmed by the qRT-PCR results, implying the reliability and accuracy of the RNA-Seq analysis (Fig. 6).

Discussion
Fish exposed to thermal/cold conditions will show some signs of stress, which results in the depression of immune responses 22 , reproduction 23 , energy metabolism 7 and growth 2 . In November, the surface temperature of the northwestern Pacific Ocean was unstable, and the temperature difference in one photoperiod varied from 7 °C to 20 °C. Under this environment, the black rockfish may experience serious acute temperature stress, which will cause heat shock, disease, and metabolism and reproduction problems. However, studies investigating the molecular mechanisms under temperature stress in black rockfish are still lacking. As studies have shown that the liver is one of the most important organs for metabolism adjustments in the process of stress adaptations 24 , in the present study, we conducted an RNA-Seq analysis on liver samples to reveal the molecular mechanisms underlying the response to temperature stress in black rockfish.  A total of 250,326 transcripts were generated with 66,596 (30.7%) transcripts yielding the Nr databases match, which greatly enriched the transcriptome data of black rockfish. This study not only identified potentially differentially expressed transcripts under acute thermal/cold conditions but also identified many new annotated gene sequences in black rockfish.
To maintain homeostasis under acute stress, energy supply and immune response pathways are activated, along with the activation of material synthesis, metabolic activity and signal pathways. In this present study, a total 584 annotated transcripts were identified in the black rockfish liver during the three comparisons (HT, LT and NT) in response to temperature stress. These differentially expressed genes were enriched and categorized based on a GO annotation, KEGG enrichment analysis and manual literature search, and several key genes or pathways likely involved in responses/adaptions to temperature stress were highlighted, as discussed below.
Candidate genes or pathways involved in the heat stress response. In this study, 8 differently expressed genes were identified by the HT vs LT vs NT comparison (Fig. 4): 1 gene related to signal transduction (early growth response protein 1), 1 gene related to molecule transport (bile salt export pump, abcb11), 1 gene related to protein folding (hsp70a), 1 gene related to immune responses (rtp3), 2 genes involved in metabolism (1,25-dihydroxyvitamin d(3) 24-hydroxylase, apoa4), and 2 genes involved in other functions (unchartered gene, transcription factor jun-b-like). Therein, HSP70 is a charter stress response gene and has been mentioned along with Apoa4 in a previous study. Heat shock stress is considered to be a well-known and studied stressor. In a study on grass carp (Ctenopharyngodon idellus), HSP70 gene expression was found to be up-regulated in spleens under high temperature stress 22 . Early growth response protein 1 (EGR1) indicates that the DNA methylation status of the promoter under stress plays a crucial role in the consolidation of immobility behavior 25 . Transcription factor jun B is a part of the inducible transcription factor complex AP-1, which is quickly activated during gravity alterations and regulates the formation of primary osteoblasts 26 . Bile salt export pump (ABCB11) functions in bile acid transport and is a susceptive factor in hepatocytes injury 27 , and, similar to the results observed here, it was reduced after heat stress during a previous study on rats 28 . Receptor transporting protein 3 (RTP3) was found to be associated with virus infection in Asian seabass 29 .
Among the 9 KEGG pathways enriched in both the HT vs LT group and the HT vs NT group, 4 KEGG pathways (influenza A, legionellosis, lnflammatory bowel disease (IBD) and measles) were related to immune and infectious diseases, and 4 KEGG pathways (NOD-like receptor signaling pathway, osteoclast differentiation, plant-pathogen interaction inflammatory bowel disease (IBD) and antigen processing and presentation) were related to the immune system. Considering that heat stress has a negative effect on the inner immune system 30,31 , it is no wonder that infection-related pathways (e.g., influenza A) and immune system response pathways (e.g., NOD-like receptor signaling pathway) are activated after stress treatments.
Influenza A viruses are the agents for a disease that can lead to high morbidity 32 (Fig. 3), and legionellosis is a disease caused by Legionella cell infection 33 . In a study on mice under chronic heat stress, the inner immune system of mice was affected and infected with the influenza virus 34 . The immune system is reported to be affected by thermal stress 30,31 , so we can infer that (1) the intracellular immune system may be damaged by thermal stress, leading to infection by some pathogens, or (2) intracellular or intercellular compounds (such as protein and deoxyribonucleic acid) may be damage, activating protein folding and degradation progresses 35 . Five DEGs (hsc70, hsp70a, hsc71, il-1β and ikkalpha) were enriched in both influenza A and legionellosis pathways, which are closely related to the immune system. Among the 5 DEGs, hsc70, hsp70a, and hsc71 are three typical stress-related genes. In Atlantic salmon (Salmo salar) and brook charr (Salvelinus fontinalis), a similar response, that hsp70 protein levels increased under both thermal and cold conditions, was observed 36 . Furthermore, triploids of these two species have the same hsp70 level tendency, with a relatively low concentration compared with the diploids 36   interesting that the heat shock cognate 70-2 (hsc70) was enriched in 5 KEGG pathways (corrected p-value < 0.05): influenza a (ko05164); legionellosis (ko05134); estrogen signaling pathway (ko04915); measles (ko05162); antigen processing and presentation (ko04612), indicate that hsc70 plays a comprehensive role in the acute thermal stress response, a result that has been shown in different species 4,37,38 . Similar to the present results, other studies have reported that hsc71 is up-regulated after hypoxia stress 39 . Interleukin 1 beta (IL-1β) is an evolutionarily conserved molecule originally identified in the immune system, and it plays a critical role in the activation of immune cells 40 . Notable, a study by Tort L showed that the fish immune response is activated under acute stress, but is suppressed under chronic stress 41 . In addition, a study on the heat shock responses of rats suggest that IL-1b plays a major role in heat-induced liver damage, and plays an important role in hepatocyte apoptosis in heat-induced liver injury 42 . In this work, the potentially differentially expressed genes with critical roles in immune responses were functionally annotated (Table 3). Interleukin-1 beta (IL-1β) showed an up-regulation trend in the HT treatment, similar to the results reported by a study on the Chinese brown frog (Rana dybowskii) 40 , a vertebrate, as well as results reported by a study on the skeletal muscles of Sebastes schlegelii 43 . In addition, in this study, ikkalpha was up-regulated after acute heat stress. However, in some other heat stress studies, the ikkalpha protein was found to be depleted and phosphorylated in male Sprague-Dawley rats 44 or coprecipitated with Hsp90 45 . Thermal stress can result in serious stress-associated inflammatory and metabolic changes 46 . The main function of the NOD-like receptor signaling pathway (Fig. 4) is inflammasome activation 47 , and osteoclast differentiation has been reported to be associated with the immune system 48 . In this study, black rockfish were under acute thermal stress, and pathogen infections may activate the organismal immune system, causing serious pathway activation. It is well known that ERs participate in the transcription complex with a number of chaperones and cofactors, including HSPs 49 . ER-binding HSP90 is accessible for hormone binding; furthermore, hormone binding promotes a receptor isolated from HSP90, converting it into a DNA-binding state 50 . The sugt1 protein has been shown to be a binding partner of heat shock proteins, and has been found to increase after heat stress 51 , and sugt1, along with hsp90, is found to be essential in both mammalian and plant innate immune responses 52 . In the 4 pathways related to the immune system mentioned above, c-jun was up-regulated in the osteoclast differentiation and IBD pathways, and a similar result was found in mice skeletal muscles after heat stress 43 . It has been suggested that c-jun may participate in signal transcription to induce an early stress-induced immune response 43 .
Four hsp genes and 1 pin1 gene related to protein folding were enriched in the estrogen signaling pathway in the liver of black rockfish after acute thermal stress (Fig. 5). A similar hsp70/90 up-regulating response was observed in rainbow trout 53 , which suggests that hsp90 is necessary for vitellogenin induction which is the production of the estrogen signaling pathway.  Metabolism. Temperature changes may influence aspects of metabolism especially in the oxidant reduction process 6 , and responses to stress are energy-costing processes 7 . In our results, some differently expressed genes involved in glycogen synthesis, fatty acid synthesis and oxidant reduction were overexpressed in the liver.
Malate dehydrogenase (MDH) is one key enzyme in the conversion of malate and oxaloacetate by the NAD/NADH system. It is well known in kinetic studies that NAD/NADH is the first compinent in the reaction of malate to oxaloacetate 54 . Under thermal and cold stress in Sebastes schlegelii, MDH and NADH-related genes (MT-ND3, MT-ND4 and MT-ND5) were down-regulated in the HT vs LT group. The same result of heat-stress-induced repression in genes encoding enzymes (MT-ND1, MT-ND2 and MT-ND6) was revealed in channel catfish (Ictalurus punctatus) 7 . Cytochrome P450 (CYP) is a superfamily containing a series of genes encoding P450 enzymes and are found in all aerobic eukaryotes and other vertebrates 55,56 . In a study on the cytochrome p450 metabolic enzymes in cows under heat stress conditions, the relative abundances of CYP2C and CYP3A were found to be decreased 57 . A study on Symbiodinium under both rapid and gradual thermal stress revealed that up-regulation occurred under gradual heat stress before the maximum temperature was reached, and down-regulation occurred under rapid stress and gradual stress after the maximum temperature was reached 58 . In a study on Sebastes schlegelii, MDH and NADH-related genes (MT-ND3, MT-ND4 and MT-ND5) were down-regulated under acute thermal stress compared with the cold stress group and the control group, which are both involved in the tricarboxylic acid cycle, a crucial pathway in oxidative metabolism. Importantly, cyp3a4 and the associated linoleic acid metabolism pathway (ko00591) were significantly changed under thermal stress, which is similar with the results of a study investigating mitochondrial functions following hypoxia 59 . Under hypoxia caused by the thermal environment, the anaerobic metabolism level will rise, while oxidative metabolism will be repressed, which results in the down-regulation of the oxidative metabolism enzyme mentioned above. Similar effects on anaerobic metabolism caused by thermal inducement have been observed in other fish species 7,60 . Lactate dehydrogenase (ldh) and cytochrome c (cyt) were observed to be increased under acute warm conditions in a study on rainbow trout (Oncorhynchus mykiss) 61 , which agrees with the results in LT treatment of the present study, suggesting energy consumption and functional impairment in mitochondria. In addition, some other metabolic-related genes, such as glycine dehydrogenase (gldc), insulin-induced gene 1 protein (insig1), thioredoxin reductase 1 (txnrd1) and apolipoprotein A-IV (apoa4), all showed significant changes in this study, especially under thermal stress 62,63 . Protein folding. Temperature affects protein synthesis, modification and degradation at the cellular level because proteins are denatured or misfolded and then become cytotoxic by forming aggregates 5,7,64 . With increases in the number of damaged proteins, the regulation of the repair and degradation of denatured proteins subsequently activates to maintain homeostasis in the cell 65 , a process in which some chaperone proteins are involved.
Heat shock proteins (HSPs), also known as stress proteins, are among the molecular chaperones that play a fundamental role in the regulation of normal protein synthesis and produced in all cellular organisms exposed to stress 66 . A study on in blue-green damselfish (Chromis viridis) observed that HSP70 and HSP60 were both elevated in response to a temperature of 32 °C 65 . The present study observed enrichment of HSP70, HSP90, and HSP40 family and other heat shock protein genes in Sebastes schlegelii, with a significant elevation of all these HSPs observed under acute heat stress (Table 3). HSP40 plays a role in the regulation of HSP70 activity by interacting with both the HSP70 and J domains 7 . HSP90 is essential in the folding and assembly of cellular proteins and is involved in the regulation of kinetic partitioning among folding, translocation and aggregation in the cell 66 , especially under damages caused by thermal and cold stress conditions. Protein modification, by folding degradation represents a series of complex pathways involving different molecules. Ubiquitin in cells acts as a covalent modifier of proteins in functionalization and degradation, which is dependent on ubiquitin ligase. E3 ubiquitin proteins are the final enzymes in the ubiquitin-proteasome pathway, regulating protein degradation, cell growth and apoptosis in response to environmental accommodation 67 . In addition, stress-induced phosphoprotein 1 (stip1) is also known as an HSP70/HSP90 organizing protein, expressed in the heat shock response 68 .

Signal transduction.
Responses and accommodations to different stresses involve a series of comprehensive and complex pathways. G protein-coupled receptor 155 (GPR155) belongs to the seven-transmembrane domain of the GPCRs superfamily 69 . The ligands for the GPCRs have varied ions, amines, proteins and lipids, which may be caused by stress and accommodation 70 . The CREB (cAMP response element binding) protein is a cellular transcription facto that responds to different physiological signals, including stresses 9 . In this study, some differentially expressed genes potentially involved in signals transduction were found in the HT vs LT group, such as G protein-coupled receptor 155 (GPR155), MAP kinase-interacting serine/threonine-protein kinase 2 (MNK2), methionine tRNA ligase and cyclic-AMP response element-binding protein 2 (CREB). They may have important functions in regulation of signaling to activate responses against harm caused by thermal conditions. Immune response. Different from mammals or birds, fish are ectothermic, with immune systems exposed to changes in the external temperature 71 . In addition, teleost have a complete vertebrate immune system similar to that of mammals 72 . Previous studies have focused on immune responses within different temperatures ranges in different fish species. Complement C3 protein gene expression increased in the orange spotted grouper (Epinephelus coioides) liver under temperature stress, and C3 may play a critical role in immune mechanisms 73 . Some other genes were found to be potentially involved in the S.schlegelii heat stress response, such as the C-X-C motif chemokine 11 (CXCL), which is an interferon-induced inflammatory chemokine expressed by leukocytes, fibroblasts and endothelial cells 74 , latexin (Lxn) and complement C3. Further immune response mechanisms will studies in more detail in the future.
In conclusion, the results of this study demonstrate that the the acute thermal conditions and stress significantly affect black rockfish (S. schlegelii). A total of 584 DEGs were obtained in response to acute thermal (27 °C) and cold (5 °C) stress exposure, such as hsps, mdh, cyp2c. These stress-regulated genes are associated with metabolism, protein folding, immune response, cell proliferation/apoptosis, membrane, molecule transport, regulation of transcription and others categories, which enables the understanding of molecular mechanism in response to temperature stress for aquatic species.

Materials and Methods
Ethics statement. All procedures involved in handling and treatment of fish during this study were approved by Animal Research and Ethics Committees of Ocean University of China prior to the initiation of the study. The field studies did not involve endangered or protected species. All experiments were performed in accordance with relevant guidelines and regulations. Animals. 40 male adults of black rockfish cultured by cages were obtained in November from northern Yellow Sea, Shandong province, China. The natural seawater temperature was 16 °C (±0.5 °C). Following capture, fish were acclimatized at a density of 10 individuals per tank (diameter 1 m, height 1.5 m) under laboratory conditions for two days without feeding. Water temperature, dissolved oxygen and salinity were maintained at 16 °C (±0.7 °C), 7.22 mg/L (±0.59 mg/L) and 30 ppm, respectively. Temperature challenge and fish sampling. After acclimation, a total of 30 fishes were randomly divided into 3 groups: low temperature group (LT, n = 10), control group (natural temperature, NT, n = 10) and high temperature group (HT, n = 10). The temperatures of the above three groups were set at 5 °C (±0.5 °C), 16 °C (±0.5 °C) and 27 °C (±0.5 °C), respectively. Three water tanks were filled by fresh seawater which was heated using heating rod or cooled down by refrigerator before treating.
Fish were transferred to the three water tanks directly by groups. After 12 h treating, 10 individuals per tank were sampled under 200 mg /L tricaine methanesulfonate (MS-222) anesthesia. Liver samples were collected from all individuals in each treatment, which were frozen immediately by liquid nitrogen and stored at −80 °C for RNA extraction.
RNA extraction, library construction and transcriptome sequencing. To reduce the variation among individuals, 3 liver tissue samples/ treatment were mixed for RNA extraction. Total RNA was extracted from freshly thawed liver samples using TRIzol ® reagent (Invitrogen, USA) and treated with TURBO DNA-free ™ kit (Invitrogen) to remove genomic DNA. The concentration and quality of the total RNA were assessed by Agilent 2100 Bioanalyzer system (Agilent Technologies, USA). Further, equal volume of RNA from 3 mixed samples/ treatment were pooled together in order to mask the difference among repetitions. 6 sequencing libraries totally were generated using NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following manufacturer's instructions and index codes were added to attribute sequences to each sample. Samples were sequenced on an Illumina Hiseq. 4000 platform and 150 bp paired-end reads were generated. Raw sequences were deposited in the Short Read Archive of the National Center for Biotechnology Information (NCBI) with accession numbers of SRR4409372 (NT), SRR4409389 (LT) and SRR4409390 (HT).
Quality control and De novo assembly of sequencing reads. Initially, reads with adapter, reads containing more than 0.1% poly-N and low quality reads were trimmed to generate high quality clean data. Then de novo assembly was performed on liver clean reads using the Trinity assembly software suite 75 . Trinity's assembly pipeline consists of three consecutive modules: Inchworm, Chrysalis, and Butterfly. All overlapping k-mers (k-mer = 25) were extracted from clean reads. Inchworm then examined each unique k-mer and generated transcript contigs using a greedy extension based on (k-1)-mer overlaps. Chrysalis clusters related Inchworm contigs into components, which were encoded by building a de Bruijn graph for each cluster. This clusters together regions that have likely originated from alternatively spliced transcripts or closely related gene families. Finally the Butterfly module processed the individual graphs in parallel, generating final transcripts 76 . Annotations of transcripts and pathways. Transcripts (both contigs and singletons) were annotated by performing BLASTx searches 77 using NCBI non-redundant (Nr), NCBI nucleotide sequences (Nt) and Swiss-Prot databases with a cutoff "e-value" of <1e-5. Domain-based comparisons with Pfam (Protein family) and KOG (a eukaryote-specific version of the Clusters of eukaryotic Ortholog Groups) databases were performed by RPS-BLAST tool from locally installed NCBI BLAST + v2.2.28 and HMMER 3.0 program, respectively. Annotated transcripts were analyzed to gene ontology (GO) classification with the aid of Blast2Go program 78 . These gene terms were then enriched on the three GO categories (Biological Process, Cellular Component and Molecular Function at level 2) using the GOseq R package 79 . Kyoto Encyclopedia of Genes and Genomes (KEGG), which is a database of biological systems, maps were retrieved by online KEGG Automatic Annotation Server for the overview of metabolic pathway analysis 17,19,20,80 . Differential gene expression analysis. The reads of each library were separately mapped to the de novo assembled transcripts with the aid of bowtie 2 program with no mismatch 81 . Count numbers of mapped reads and FPKM (expected number of Fragment Per Kilobase of transcript sequence per Millions base pairs sequenced) were retrieved and normalized by RSEM V1.2.15 82 . Differential expression statistical analysis of three treatment (NT, LT and HT) was conducted by the DEGSeq R package 83 with a cutoff "q-value" of 0.05 and |log2(fold change)| > 1. Transcripts with absolute fold change values over 2.0 were marked as significantly differential expressed genes.
Experimental validation by quantitative real-time PCR. 11 differentially expressed genes were randomly selected for validation using quantitative real-time PCR (qRT-PCR) with gene specific primers designed using Primer 5 software (Premier Biosoft International) to validate our Illumina sequencing data. Primers were listed in Supplement 6. Samples were generated from NT, LT and HT groups in the preceding experiment. The first strand cDNA was synthesized by using M-MLV Reverse Transcription Kit (Promega, USA) from 1 μg of RNA. All the cDNA products were diluted to 500 ng/μl. The 20 μl qRT-PCR reaction mixture consisted of 2 μl cDNA template, 0.4 μl of both primer, 10 μl of KAPA SYBR ® FAST qPCR Master Mix (2×), 0.4 μl of ROX and 6.8 μl of RNAase-free water. PCR amplification was performed as that incubated in a 96-well optical plate at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 58 °C for 30 s, and a final extension at 72 °C for 2 min. qRT-PCR was performed using the StepOne Plus Real-Time PCR system (Applied Biosystems) and 2-ΔΔCT method was used to analysis the expression level of genes. 18S ribosomal RNA (18S) and were used as the reference gene for qRT-PCR normalization.