Transcriptomic analysis of the hepatic response to stress in the red cusk-eel (Genypterus chilensis): Insights into lipid metabolism, oxidative stress and liver steatosis

Teleosts exhibit a broad divergence in their adaptive response to stress, depending on the magnitude, duration, and frequency of stressors and the species receiving the stimulus. We have previously reported that the red cusk-eel (Genypterus chilensis), an important marine farmed fish, shows a physiological response to stress that results in increased skeletal muscle atrophy mediated by over-expression of components of the ubiquitin proteasome and autophagy-lysosomal systems. To better understand the systemic effects of stress on the red cusk-eel metabolism, the present study assessed the transcriptomic hepatic response to repetitive handling-stress. Using high-throughput RNA-seq, 259 up-regulated transcripts were found, mostly associated with angiogenesis, gluconeogenesis, and triacylglyceride catabolism. Conversely, 293 transcripts were down-regulated, associated to cholesterol biosynthesis, PPARα signaling, fatty acid biosynthesis, and glycolysis. This gene signature was concordant with hepatic metabolite levels and hepatic oxidative damage. Moreover, the increased plasmatic levels of AST (aspartate aminotransferase), ALT (alanine aminotransferase) and AP (alkaline phosphatase), as well as liver histology suggest stress-induced liver steatosis. This study offers an integrative molecular and biochemical analysis of the hepatic response to handling-stress, and reveals unknown aspects of lipid metabolism in a non-model teleost.

Introduction Stress generated during intensive fish farming is a major factor negatively affecting animal growth and health, and results in a decline in finfish aquaculture production [1,2]. Among the main types of stressors are those induced by physical stimuli, such as transport, capture, or handling [3,4], which are unavoidable in intensive aquaculture [5]. Fish possess a set of physiological strategies, allowing them to respond to stressors through adaptive neuroendocrine adjustments, collectively termed the stress response [6]. In fish, this stress response is mediated by the hypothalamic-pituitary-interrenal (HPI) axis, which is responsible for promoting the cortisol synthesis and secretion through interrenal cells into the bloodstream [7].
The liver is one of the most important target organs in the adaptive response to stress, and is a fundamental metabolic tissue for energy substrate administration [8], because hepatic tissue can synthesize de novo glucose for non-hepatic tissues during periods of stress [9]. In addition, the liver is a key tissue involved in somatic growth regulation, the immune response, and detoxification [10][11][12][13]. Recent advances in transcriptomic technologies now allow a comprehensive understanding of how physical disturbances modulate metabolic adaptations to stress [14,15]. While early studies were conducted on a limited number of species using microarray technology [16][17][18], it is now possible, thanks to progress in high-throughput sequencing technologies, to study the adaptive function of liver tissue in a range of commercially important species, including channel catfish (Ictalurus punctatus) [19], rainbow trout (Oncorhynchus mykiss) [20], and large yellow croaker (Larimichthys crocea) [21], among others. Such studies demonstrate that each species has different adaptive mechanisms and molecular responses to stress, revealing particular details in their species-specific stress tolerances, meaning that culturing a new species requires an integrative and detailed knowledge of tissue-specific stress responses.
The red cusk-eel (Genypterus chilensis) is a new economically important marine species for the Chilean aquaculture industry [22], but low handling-stress tolerance results in high mortality rates during juvenile stages [23]. Handling-stress facilitates red cusk-eel skeletal muscle atrophy through the coordinated expression of components of the ubiquitin-proteasome and autophagy-lysosome systems, revealing particularities in the compensatory response to stress [24]. Therefore, in the present study, we evaluate the hepatic response to handling-stress in the red cusk-eel using RNA-seq analysis, using a recently reported hybrid reference transcriptome [25]. This transcriptomic analysis was complemented by a quantification of hepatic metabolites, plasmatic enzyme activities, liver histology and oxidative stress markers. Our results strongly suggest that handling-stress has detrimental effects in this species, modulating the expression of genes associated with hepatic lipid metabolism and inducing liver steatosis.

Animals and experimental design
Juvenile red cusk-eels (Genypterus chilensis) with an average weight of 900 ± 50 g and length of 55 ± 5 cm were collected from the Centro de Investigación Marina de Quintay (CIMARQ) (33˚13 0 S 71˚38 0 W, Valparaíso Region, Chile). Eight fish were kept under natural temperature and light:dark photoperiod conditions (13˚C ± 1˚C and L:D 12:12) for the spring season. The animals were then split into control and stressed groups, placed in separate 90 L tanks, and acclimated for two weeks before the handling-stress protocol. Fish were fed once daily with 6-mm commercial pellets, containing 60% protein, 6% lipids, 11% carbohydrates, 15% ashes, and 8% humidity (Skretting, Puerto Montt, Chile). The stressed group was then subjected to a standardized handling-stress protocol consisting of netting and chasing the fish for five minutes daily for five days. The management bias from adjacent tanks was prevented through the completely blocked by black covers on all tanks. Six hours after the final handling stimulation, the four stressed and four control fish were captured and anaesthetized (3-aminobenzoic acid ethyl ester, 100 mg/L -1 ). Blood was collected from the caudal vein with 1 ml heparinized (10 mg/mL) syringes and centrifuged at 5000 xg for 10 min at 4˚C to obtain plasma, frozen in liquid nitrogen, and stored at -80˚C. The fish were then euthanized through an overdose of anesthetic (3-aminobenzoic acid ethyl ester, 300 mg/L). Livers were collected and immediately frozen in liquid nitrogen and stored at -80˚C.

AST, ALT and AP measurements
The plasmatic activity of AST (aspartate aminotransferase), ALT (alanine aminotransferase) and AP (alkaline phosphatase) were determined using commercially available kits from Valtek (Santiago, Chile). The enzymatic activity of AST, ALT, and AP are determined, respectively, by the production of colorimetric products from the generation of glutamate (colorimetric product: 450 nm), pyruvate (535 nm), and p-nitrophenol (405 nm).

Hepatic metabolites
The hepatic metabolite levels for glycogen, triglyceride, non-esterified free fatty acid (NEFA) and total cholesterol were determined using commercially available kits from Cell Biolabs (CA, USA). Ultrasonic disruption in PBS containing 0.5% Triton X-100 was used to homogenize 200 mg of liver tissue, which was centrifuged at 5000 xg for 10 minutes at 4˚C. The supernatant was collected, stored on ice, and used for glycogen and triglyceride quantification using the Glycogen Assay Kit (Cell Biolabs), and Serum Triglyceride Quantification Kit (Cell Biolabs), respectively. For cholesterol and free fatty acid quantification [26], 10 mg of liver tissue was extracted with 200 μL of a chloroform:isopropanol:NP-40 (7:11:0.1) mixture in a microhomogenizer. The extract was centrifuged at 15,000 xg for 10 minutes, and the organic phase was transferred to a new tube and air dried at 50ºC to remove the chloroform. The dried lipids were dissolved and homogenized in 200 μL of 1X Assay Diluent via vortexing. Cholesterol and NEFA were quantified using Total Cholesterol Assay Kit (Cell Biolabs) and Free Fatty Acid Assay Kit (Cell Biolabs), respectively.
Hepatic oxidative stress assays DNA oxidative damage, protein carbonylation, and lipid peroxidation were determined using commercially available kits from Cell Biolabs (CA, USA). For oxidative DNA damage assessment, genomic DNA (gDNA) was extracted from 25 mg of liver using the Isolate II Genomic DNA Kit (Bioline, MA, USA), and dissolved in TE buffer. DNA was quantified by spectrophotometry using NanoDrop with the Epoch Multi-Volume Spectrophotometer System (BioTek, VT, USA). Assessment of gDNA quality and visualization were performed by electrophoresis on a 1% agarose gel in TAE 1X containing ethidium bromide. In order to determine apurinic/ apyrimidinic (AP) sites, gDNA was mixed with the aldehyde reactive prove (ARP) reaction using the OxiSelect Oxidative DNA Damage Quantification Kit (Cell Biolabs, CA, USA), reading absorbance at 450 nm. For assessment of protein carbonylation, total proteins were extracted from 100 mg of liver in 1 ml of lysis buffer containing 50 mM Tris-HCl pH 7.4, 150 mM NaCl, 1 mM EDTA, and 1% NP-40, and solubilized at 4˚C after 12,000 xg centrifugation. Protein concentration was determined using a Pierce1 BCA Protein Assay Kit (Thermo Scientific, IL, USA). Then, protein carbonyl content was quantified by the OxiSelect Protein Carbonyl Spectrophotometric Assay (Cell Biolabs, CA, USA). Lipid peroxidation assessment was quantified by the formation of thydroxynonenal (HNE) protein adducts using the OxiSelect HNE Adduct Competitive ELISA Kit (Cell Biolabs, CA, USA).

Liver histology
The livers of the red cusk-eel were dissected immediately after euthanization according to the protocol described previously. The livers were fixed in 4% paraformaldehyde, embedded in paraffin, and sectioned at a 5 μM thickness of three slices per liver. The samples were stained with hematoxylin/eosin, observed in a Olympus BX-61 microscope and photographed with a Leica DF300 camera.

Liver RNA sequencing
Total RNA was extracted using the RNeasy Mini Kit (Qiagen, TX, USA). RNA was quantified spectrophotometrically using NanoDrop with the Epoch Multi-Volume Spectrophotometer System (BioTek, VT, USA). Total RNA isolated from the liver was treated with DNase I to remove genomic DNA. RNA concentration was measured using a Qubit1 2.0 Fluorometer (Life Technology, Carlsbad, CA, USA), and RNA integrity was determined using the Fragment Analyzer™ Automated CE System (Analytical Advanced Technologies, Ames, IA, USA). Equal quantities of total RNA from four fish were pooled by condition, and the total RNA was used to prepare mRNA libraries. Complementary DNA (cDNA) libraries were constructed using Tru-Seq RNA Sample Preparation kit v2 (Illumina1, USA). Libraries were sequenced (2 x 250bp) with the MiSeq (Illumina1) platform at the Centro de Biotecnología Vegetal (UNAB, Chile).

Analysis and representation of differentially expressed transcripts
Raw sequencing reads were trimmed by removing adaptor sequences, low quality sequences (quality scores less than 10), and sequences with lengths less than 30 bp. To identify differentially expressed transcripts, the reads from control and stressed conditions were mapped to the G. chilensis reference transcriptome [25], using CLC Genomics Workbench, v.7.0.3 (http:// www.clcbio.com), with the following parameters: mismatches = 2, minimum fraction length = 0.9, minimum fraction similarity = 0.8, and maximum hits per read = 5. Gene expressions were based on reads per kilobase of exon model per million mapped read (RPKM) values. Transcripts with absolute fold-change values > 2.0 and FDR corrected P value < 0.05 were included in the GO and KEGG enrichment analyses.

RNA-seq validation by real-time qPCR
All quantitative real-time polymerase chain reaction (qPCR) assays were carried out complying with the MIQE guidelines [28]. Total RNA was extracted from liver tissue using the RNeasy Mini Kit (Qiagen, TX, USA). RNA was quantified by spectrophotometry using Nano-Drop with the Epoch Multi-Volume Spectrophotometer System (BioTek, VT, USA). Only RNA with an A260/280 ratio between 1.9 and 2.1 were used for cDNA synthesis. Residual genomic DNA was removed using the genomic DNA wipeout buffer included in the Quantitect Reverse Transcription Kit (Qiagen, TX, USA). Subsequently, 1 μg of RNA was reverse transcribed into cDNA for 30 min at 42˚C.
The qPCR assessments were performed using a Stratagene MX3000P qPCR system (Stratagene, La Jolla, CA, USA). Each reaction mixture contained 7.5 μl of 2× Brilliant1 II SYBR1 master mix (Stratagene, La Jolla, CA, USA), 6 μl of cDNA (40-fold diluted) and 250 nM of each primer in 20 μl final volume. The list of primers and amplification efficiencies are described in S1 Table. Amplifications were performed in triplicate following thermal cycling conditions of: initial activation of 10 minutes at 95˚C, followed by 40 cycles of 30s of denaturation at 95˚C, 30s of annealing at 54-60˚C, and 30s of elongation at 72˚C. In order to confirm the presence of a single PCR product, dissociation curve analysis of the PCR products was performed and evaluated by electrophoresis on a 2% agarose gel. With the purpose of calculating the assay's efficiency, 2-fold dilution series were created from a cDNA pool. Efficiency values were estimated from the slope of the efficiency curve: E = 10(-1/slope) -1 . The reference gene, 40S ribosomal protein S30 (fau), was used for gene expression normalization, and control reactions included a no-template control and a no-reverse-transcriptase control. QGene was used for analysis of gene expression [29].

Statistical analysis
Data are expressed as the mean ± SEM. Differences in means between groups were determined using one-way ANOVAs, followed by a Bonferroni post hoc test. Data were accepted as significant at a value of P < 0.05. Correlations between RNA-seq and qPCR data were assessed through multiple linear regression, using coefficients of determination (R 2 ) and p-values. All statistical analyses were performed using GraphPad Prism, v.5.00 (GraphPad Software, CA, USA).

Assessment of hepatic metabolites
We previously reported that plasmatic cortisol and glucose levels significantly increased after 5 days of daily handling-stress [24]. Considering how the liver plays a central role in the metabolic homeostasis of carbohydrates and lipids, we measured the hepatic metabolite levels of glycogen, triglyceride, cholesterol and non-esterified fatty acids (NEFA). We found decreased glycogen (Fig 1A), triglyceride ( Fig 1B) and cholesterol levels ( Fig 1C) in comparison to control conditions. Conversely, the levels of NEFA were higher than in the control condition ( Fig  1D). These results indicate a dynamic response to stress-mediated carbohydrate and lipid metabolism.

Hepatic transcriptomic responses of red cusk-eel to handling-stress
To understand how stress modulates the hepatic stress response in a comprehensive manner, RNA-seq analyses were performed on liver samples from control and stressed groups including cDNA library replicates of each condition. We obtained a total of 12,950,452 paired-end reads for the control and a total of 13,669,452 paired-end reads for the stressed group. Raw reads have been deposited in the NCBI Sequence Read Archive under the study accession number SRX1056881. After sequence trimming for adapter subtraction and filtering of low quality base pairs, the data sets were reduced to 11,912,668 and 12,270,950 high-quality reads, respectively. Differentially expressed transcripts (DETs) were estimated by FPKM, mapping obtained reads to the reported reference transcriptome for G. chilensis [25], resulting in a mapping of roughly 98.5% of the reads. A total of 552 DETs showed significant differential expression between both conditions. Of these, there were 259 up-regulated transcripts and 293 down-regulated transcripts in stressed fish. A complete list of the differentially expressed transcripts is included in S2 Table. GO cluster enrichment and pathway analysis The up-and down-regulated genes were analyzed through the DAVID database and clustered based on GO terms. Table 1 summarizes the top four clusters of functional categories with Enrichment Score > 4 up-regulated under stress. The cluster with the highest enrichment scores (ES 6.69) includes biological processes like vasculature development, blood vessel development, blood vessel morphogenesis, and angiogenesis. Other enriched clusters up-regulated by stress included biological processes like cell migration, molecular function like heparin binding and biological process like response to hormone stimulus. Table 2 summarizes the top four clusters of functional categories with Enrichment Score > 4 down-regulated under stress. The cluster with the highest enrichment scores (ES 17.02) includes biological processes like steroid metabolic process, sterol metabolic process, sterol biosynthetic process, cholesterol metabolic process, cholesterol biosynthetic process, steroid biosynthetic process, lipid biosynthetic process. Other enriched clusters down-regulated by stress included cellular component like endoplasmic reticulum, molecular function like monocarboxylic acid binding and cellular component like microsome.

Hepatotoxicity induced by handling-stress
To evaluate the hepatotoxicity induced by this stress, we measured cellular and plasmatic markers of liver damage. Handling-stress significantly increased the rate of lipid peroxidation (Fig 3A), protein carbonylation (Fig 3B), and DNA oxidative damage (Fig 3C) (3.4-fold, 2.8-fold, and 3.3-fold, respectively) in comparison to control conditions. Similarly, the plasmatic concentrations of alanine aminotransferase (Fig 4A), aspartate aminotransferasa ( Fig  4B) and alkaline phosphatase (Fig 4C) were significantly raised in stressed animals (1.3-fold, 3.8-fold, and 1.8-fold, respectively), in comparison to control conditions. In addition, histological liver analysis reveals that hepatocytes from control conditions have the typical polygonal cell shape and clear cytoplasm. However hepatocytes from stressed conditions start to lose their polygonal shape and have a condensed cytoplasm characteristics of liver fibrosis (S7 Fig).

Discussion
Our results are the first to describe the molecular and biochemical hepatic response to stress in the marine teleost, the red cusk-eel. This species' molecular response to stress has been documented in terms of the effects of cortisol-induced stress in skeletal muscle atrophy [24]. This work provides new evidence of the susceptibility of this non-domesticated organism and reinforces the idea that each species has specificities in their stress tolerance. Here, RNA-seq analysis evaluated the hepatic transcriptomic response to handling-stress, revealing a differential gene expression primarily associated with biological pathways, such as angiogenesis, fatty acid biosynthesis, cholesterol biosynthesis, PPARα signaling, triglyceride metabolism, gluconeogenesis and glycolysis. The results obtained by this in silico approach were validated by qPCR analysis of representative cholesterol and angiogenesis pathway genes, as well as the assessment of hepatic metabolites and markers of hepatic damage.
https://doi.org/10.1371/journal.pone.0176447.g002 phase, including phosphoenolpryruvate carboxykinase (pck1) and glucose-6-phosphatase (g6pc), as well as a decrease in the expression of pyruvate kinase (pklr) and glucokinase (gck) involved in glycolisis [32,33]. Most expressions of these genes are directly regulated by the glucocorticoid receptor [9,14]. Consistent with these observations, handling-stress in the red cusk-eel induced an up-regulation in the expression of gluconeogenic gene g6pc and a downregulation of glycolitic genes pklr and gck. These results are similar to those of juvenile carp (Cyprinus carpio) and Mozambique tilapia (Oreochromis mossambicus), for which the administration of a high dose of cortisol significantly increased glucose serum levels and glucose-6 phosphatase activity, revealing a major role for cortisol in gluconeogenesis [34,35]. Similar to our observations, chronic cortisol administration in gilthead sea bream (Sparus aurata), induces a down-regulation of enolase 1 (eno1) [18], a glycolitic gene. In addition, in Senegalese sole (Solea senegalensis), repeated handling-stress induces a decrease of enolase 1 protein expression [36]. Therefore, we hypothesize that the increase in plasma glucose levels seen in the red cusk-eel after stress may be associated with enhanced muscle protein catabolism necessary for liver gluconeogenesis, even though increases in plasmatic levels of lactate during stress were not found [24]. The increase in plasmatic glucose levels can also be related to glycogen hydrolysis (glycogenolysis), evidenced by its decrease in the liver. Indeed, stress and cortisol have been demonstrated to play a fundamental role in glycogen metabolism in rainbow trout [37,38], although we did not find changes in the expression of genes associated with such catabolism.

Stress and lipid metabolism
The liver also plays a central role in lipid metabolism [39], and in the current study, handlingstress induced a general down-regulation in the cholesterol and fatty acid biosynthesis pathways. The stress-induced inhibition of cholesterol biosynthesis is in line with previous data on rainbow trout, where plasma cholesterol levels decreased after 30 days of high stocking densities [40]. Additionally, the down-regulation in fatty acid biosynthesis is concordant with observations of rainbow trout following handling-stress, in which fatty acid synthase (fasn) expression and activity decreased [37].
It is well known that under stressful conditions another metabolic fuel, ketone bodies, is generated through NEFA beta oxidation [40]. A major NEFA source is the hydrolysis of liver and adipose tissue triglycerides through lipase-catalyzed lipolysis [41]. Consistent with this, handling-stress induced an increase in the expression of liph and pnpla2 lipases as well as, a general down-regulation in the expression of genes associated with triglyceride synthesis. Lipase H (LIPH) is an evolutionary conserved phospholipase which catalyzes the production of fatty acids and lysophosphatidic acid (LPA) [42]. PNPLA2 is a member of the patatin-like phospholipase domain-containing protein family, which catalyzes the initial step in triglyceride hydrolysis [43]. Although there are reports of the over-expression of both genes in mammal tissues under pathophysiological conditions [44], no reports exist linking these genes to a stress response in teleosts. Consistent with the over-expression of liph and pnpla2 under stressful conditions, there was a decrease in triglyceride levels and an increase in NEFA hepatic levels. The dramatic increases in hepatic NEFA levels may be associated with enhanced triglyceride hydrolysis and impaired fatty acid beta oxidation, similar to decreases seen in hepatic triglyceride levels in rainbow trout caused by acute handling-stress [37]. It is therefore probable All data are represented as means ± SEM (n = 4). Significant differences between control and stressed groups are shown as * (P < 0.05) and ** (P < 0.01).

Stress and liver steatosis
Liver steatosis is defined as an accumulation of excessive triglycerides and other fats inside liver cells [46], although fish that are fed diets rich in unsaturated fatty acids are known to develop this condition [47], there are no reports relating this condition to stress. A direct relationship between glucocorticoids and fatty liver disease is described in mammals, with glucocorticoids acting as a promoter of lipid accumulation, leading to inflammation, and the potential for fibrosis [45]. This accumulation in hepatocyte diminishes mitochondrial oxidative capacity, further reducing the state of the electron transport chain complexes and stimulating peroxisomal and microsomal pathways for omega oxidation [45]. The consequent increased generation of reactive oxygen species (ROS) and reactive aldehydic derivatives causes oxidative stress [46].
We found increased expressions of cytochrome P450 2D6 (cyp2d6) and alcohol dehydrogenase class-3 (adh5), genes related to fatty acid omega oxidation and a consistent increase in lipid peroxidation, protein carbonylation, and DNA oxidation induced by stress. While there are no reports in teleost that link the fatty acid accumulation in liver to oxidative damage, it has been recently reported in tapertail anchovy (Coilia nasus) that loading stress induces lipid peroxidation, expression of apoptosis-associated markers, and increased plasmatic levels of ALT and AST [48]. While the authors associated the liver damage to TNFα-induced apoptosis, interestingly they also reported an increase of lipid metabolic genes, similar to our observations.
Recently, a study developed with a mammalian model of hepatocellular steatosis established the effects of free fatty acid and the mechanism associated to oxidative stress [49]. The authors reported that high concentration of free fatty acid induced excessive lipid accumulation and oxidative stress by down-regulating the ppara expression, similar to our observations [49]. PPARα is nuclear hormone receptors, and regulates the expression of many genes involved in fatty acid oxidation, ketogenesis, gluconeogenesis, cholesterol catabolism, among others. In fish, PPARα have been identified and characterized in several species however there are no direct reports that associate changes in its expression during stress response [50]. However, the systemic administration of the PPARα antagonist, WY-14643, has been shown to induce a decrease in NEFA plasmatic levels associated to increased fatty acid beta oxidation in turbot (Scophthalmus maximus) [51]. Therefore, we speculate that stress-impaired PPARα expression should induce an increase in NEFA overload and cytochrome P450-dependent-omega oxidation. Consistent with these observations, a direct role for glucocorticoids in the increase of TG hydrolysis, as well as in the inhibition of beta oxidation and cholesterol biosynthesis, has been demonstrated in mammals through in vitro and in vivo trials [52,53].

Conclusions
Handling-stress induced major changes in the red cusk-eel hepatic metabolic response. Under stress, 259 genes, mostly associated to angiogenesis, gluconeogenesis, and triglyceride degradation, were up-regulated. Conversely, 293 genes, associated to PPARα signaling, fatty acid biosynthesis, cholesterol biosynthesis, and glycolysis, were down-regulated. These transcriptomic analyses are supported by quantifications of liver metabolites, oxidative damage, markers of hepatic damage and liver histology. This is the first study to suggest that handling-stress is powerful inducer of liver steatosis in a teleost, providing valuable information for monitoring the culturing and growth of marine fish species under intensive rearing conditions. Supporting information S1 Table. Primer sequences for qPCR assay, amplicon size, and PCR efficiencies for genes used in the study. Sequences were derived from the red cusk-eel reference transcriptome.