Transcriptomic Analysis of Liver Tissue in Fat Greenling (Hexagrammos otakii) Exposed to Elevated Ambient Ammonia

Fat greenling (Hexagrammos otakii) is an important aquaculture fish species in northern China. Unfortunately, the increase of nitrogen in water due to residual feed and animal feces has caused considerable losses to the fat greenling breeding industry. Studies on detoxification metabolism and immune responses of fat greenling under ammonia stress have not been carried out. We used high-throughput sequencing to extract RNA from the liver of fat greenling, 48 and 96 h post-exposure, to water containing 30 mg/L NH4Cl. A total of 87,675 high-quality unigenes were obtained by transcriptome analysis, the N50 and mean length were 1,523 and 772 bp, respectively. There were 320 and 272 genes (DEGs) identified in the two sampling times (48 and 96 h), containing 132 and 137 significantly up-regulated genes, and 188 and 135 down-regulated genes, respectively. Ten DEGs were randomly selected for validation by quantitative PCR. Further annotation and analysis indicated that the DEGs were enriched in the MAPK-signaling pathway, FOXO-signaling pathway, phosphatidylinositol-signaling pathway, cytokine-cytikine receptor interaction, Hippo-signaling pathway, and neuroactive ligand-receptor interaction. These pathways were mainly related to oxidative stress and immune responses. This research provides a valuable resource to further study of detoxification metabolism and immune responses of fat greenling under ammonia stress.


INTRODUCTION
Fat greenling (Hexagrammos otakii) is one of the most important marine economic fish, largely because of its high-quality meat. In China, fat greenling are mainly distributed in the Yellow Sea, the Bohai Sea, and the East China Sea (Habib et al., 2011;Hu et al., 2017). With rapid aquaculture intensification and increased scales of production, nitrogen (especially ammonia) pollution from the aquaculture industry has become increasingly problematic. Degradation of nitrogenous organic matter, such as residuals from bait and fish excrement, increases the load of nitrogenous chemicals in surrounding waters. Increases of ammonia in aquaculture water may result in disease or death of fat greenling, which may restrict further development of the fat greenling aquaculture industry.
Ammonia nitrogen exists in two forms, namely, non-ionic ammonia (NH 3 ) and ionic ammonium (NH 4 + ); ionic ammonium is non-toxic to aquatic animals, but non-ionic ammonia may be harmful (Constable et al., 2003). Non-ionic ammonia can directly affect gills of fish and has a significant impact on enzymatic hydrolysis reactions and the stability of cell membranes. In some cases aquatic animals die, causing severe economic losses (DeFreitas et al., 2000;Koo et al., 2005). Studies also have shown that long-term exposure to NH 3 -N could affect the antioxidant enzyme activity of fish, destroy antioxidant systems, and reduce immune capacity (Hegazi et al., 2010;Yang et al., 2010Yang et al., , 2011. Under chronic ammonia stress, it also could affect osmotic pressure balance and cause death of nerve cells. As content of ammonia in tissues increases, damage to the liver and kidney systems, congestion, and hepatic comas may ensue (Felipo et al., 1994;Randall and Tsui, 2002;Eddy, 2005;Benli et al., 2008;Wilkie et al., 2011). Many studies have shown that liver plays a key role in homeostatic processes (Alvarellos et al., 2005;Alves et al., 2010;Cordeiro et al., 2012). When a fish is under stress, the liver provides energy needed by other tissues as amino acid catabolism occurs (Vijayan et al., 1990;Mommsen et al., 1999;Aluru and Vijayan, 2007), to help reduce risk of ammonia poisoning (Walsh and Mommsen, 2001;Ip et al., 2004). There has been little study of response mechanisms to ammonia stress in fat greenling. RNA-Seq is a powerful high-throughput sequencing method, which has been widely used in a variety of fish research, including minnows (Gobiocypris rarus) (Gao et al., 2015), topmouth culter (Culter alburnus) (Zhao et al., 2016), and Mandarin fish (Siniperca chuatsi) (Hu et al., 2015). To better understand the liver detoxification mechanism of fat greenling under ammonia nitrogen stress, after ammonia nitrogen exposure, we used transcriptomic profiling on fat greenling liver, targeting differentially expressed genes (DEGs) and cellular signaling pathways. We also focused analysis on the anti-oxidative system and immune responses. In doing so, we provide new data on detoxification of ammonia by fish livers.

Experiment and Sample Collection
Fat greenling (400-500 g) were obtained from Luhai Aquatic Technology Development Company of Qingdao, China. Fish were acclimated in tanks (5 m 3 ) containing aerated sand-filtered seawater (salinity 31%, pH 7.9-8.1, dissolved oxygen 6 mg/L) at 20 ± 0.5 • C for 3 weeks before experiments. During the acclimation period, half of the water was replaced twice a day in each tank, and fish were fed to satiation with special feed twice a day. Fish were fasted for 48 h before (and then during) ammonia exposure experiments.
Following studies that have previously suggested 96 h LC50 values for fat greenling was 39.35 mg/L, we performed a preliminary experiment to determine optimal ammonia concentration. Based on results from our preliminary experiment, 30 mg/L NH 4 Cl was chosen for the 48 and 96 h trials. Animals were randomly divided into the experimental group (48, 96 h) and control group (seawater only) with 3 replicates, each of which had 6 fishes. After anaesthetizing with MS-222, 3 fish livers from each replicate were removed and immediately submerged in 10 ml RNAlater (Ambion, United States) at 0 (as control), 48, and 96 h after exposure experiment. Tissues were stored at -80 • C until RNA extraction.

RNA Isolation and Illumina Sequencing
RNA was extracted from each liver sample using Trizol reagent (Invitrogen, United States) following the manufacturing protocol. RNA degradation and contamination were detected using a 1% agarose gel. RNA purity was measured using the NanoPhotometer R spectrophotometer (IMPLEN, CA, United States). RNA concentration was measured using Qubit R RNA Assay Kit in Qubit2.0 R Fluorometer (Life Technologies, CA, United States). RNA integrity was measured using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The RNA of each mixed sample is divided into two parts: one is used for transcriptome sequencing and the other is used for real-time fluorescence quantitative PCR (qRT-PCR) to validate the reproducibility of DEGs data obtained from transcriptome sequencing.
Sequencing libraries were generated using NEBNext R Ultra TM RNA Library Prep Kit for Illumina (NEB, United States) according to the manufacturer's recommendations; index codes were added to attribute sequences to each sample. Briefly, purification of mRNA from total RNA with poly-T oligo-attached magnetic beads. First, cDNA strands were synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H-). Second, cDNA strand synthesis using DNA polymerase I and RNase H. Select 150-200 bp cDNA fragments and use the AMPure XP system to purify the library fragments (Beckman Coulter, Beverly, CA, United States). The PCR products were purified by AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100 system. Then, sequencing was carried out by BioMarker (Beijing, China) using the Illumina Hiseq 2000 platform.

Transcriptome Assembly, Quality Control, and Gene Functional Annotation
The protocols of transcriptome assembly and annotation were provided by BioMarker (Beijing, China). Briefly, raw reading in fastq format is handled by an internal perl script. Clean reads were obtained by deleting reads including adapters or ploy-N, as well as low quality reads. Calculate the Q20, Q30, GC-content, and sequence duplication level of the clean data. All downstream analyses were based on high-quality cleaning data. The clean reads were spliced and assembled into transcript using Trinity software (Haas et al., 2011), and the longest transcript under each gene was extracted as a unigene (Altschul et al., 1990).

Quantification of Gene Expression Levels and Differential Expression Analysis
Gene expression levels were estimated by RSEM for each sample (Li and Dewey, 2011). The read count for each unigene comes from the mapping results and then normalized to FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) (Trapnell et al., 2010).
Differential expression was performed using the DESeq R package (1.10.1) with biological replicates (Anders and Huber, 2010). DESeq provides statistical routines for determining differential expression of digital gene expression data using a model based on the negative binomial distribution. The resulting p-value uses Benjamini and Hochberg's method to control the false discovery rate. Genes with an adjusted p < 0.05 after DESeq adjustment were assigned as differentially expressed (Benjamini and Hochberg, 1995).

GO/KEGG Enrichment Analysis for DEGs
Gene topological (GO) enrichment analysis was performed on differentially expressed genes (DEGs) using topGO R package based on Kolmogorov-Smirnov test. All DEGs of fat greenling liver were mapped to GO terms in the GO database. The study set represented the GO terms of all identified DEGs. The calculated p-value was corrected by Bonferroni, with p ≤ 0.05 as the threshold.
According to KEGG database resources, pathway enrichment analysis was performed to find significantly enriched signal transduction or metabolic pathways. KOBAS (Mao et al., 2005) software was used to detect the enrichment of differentially expressed genes in the KEGG pathway.

Quantitative Real-Time PCR (qRT-PCR)
To validate reliability of DEGs data obtained from the transcriptome sequencing in liver of fat greenling, selected 10 differently expressed genes for qRT-PCR detection. Primers for qRT-PCR were designed with Premier 5.0 (Table 1). Briefly, the same RNA samples were used for sequencing, and the cDNA was synthesized using the TransScript II All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (One-Step gDNA Removal). The qRT-PCR was carried out using Bio-Rad CFX96. All reactions were done in triplicate. The 30524 gene was used as the internal reference. Expression levels were quantitatively analyzed using the 2 − CT method (Schmittgen and Livak, 2001), and all data were given in terms of relative mRNA expression. Differences between the experimental group and the control group were used to assess changes in gene expression.

Sequencing and de novo Assembly of Fat Greenling Transcriptome Data
In total, 195,335 transcripts with an average length of 1,647bp were generated. After assembling clean reads, 87,675 unigenes were generated from 195,335 transcripts, with a N50 and mean length of 1,523 and 772 bp, respectively, were generated from 195,335 transcripts ( Table 2). There were 8,007 unigenes larger than 2,000 bp (Figure 1). The sequencing data of each treatment and control group have been stored in NCBI sequence read archive database with the accession number SAMN11155636. The major characteristics of each library, including clean reads, GC content, and Q30, are summarized in Table 3.

Functional Annotation of Unigenes
Functional annotation of unigenes was performed using the Blast comparison database. The annotated database includes COG, GO, KEGG, KOG, Pfam, Swissprot, eggNOG, and NR. A total of 28,062 unigenes were successfully matched in these databases, accounting for 32% of all unigenes (     mechanisms, post-translational modification, protein turnover, chaperones, and transcription (Figure 2). In the eggNOG database, a total of 26,230 unigenes (29.92%) were assigned to 26 descriptions. Among them, genes involved were signal transduction mechanisms, post-translational modification, protein turnover, chaperones, and intracellular trafficking, secretion, and vesicular transport (Figure 3). After GO annotation, the genes divided into three ontologies: biological process (BP), cellular component (CC), and molecular function (MF). There were 14,559 unigenes (16.6%) in the GO database which aligned with 58 functional terms. In the CC designation, these unigenes were related to the categories: cell, cell part, organelle, and membrane. In the MF designation, these unigenes mainly related to binding and catalytic activity. The BP category related to cellular process, metabolic process, singleorganism process, and biological regulation (Figure 4).

Identification and Analysis of DEGs
Analysis revealed 320 and 272 genes showed different expression in the two sampling times (48 and 96 h). Among them, the significantly up-regulated unigenes were 132 and 137, and the significantly down-regulated unigenes were 188 and 135, respectively (Figure 5). The DEGs of control group and 48 h (G0), the DEGs of control group and 96 h (G2), the DEGs of 48 and 96 h (G6) (Figure 6).

GO Functional Classification of DEGs
Go enrichment analysis of DEG showed that these genes were significantly enriched in three main functional categories (based on a value of p ≤ 0.05). GO-enrichment categories at 48 and 96 h are similar (Figure 7). In the BP category, the DEGs were significantly enriched in metabolic process, singleorganism process, and cellular process. In the CC category, the DEGs were significantly enriched in cell part and cell. In  MF category, the DEGs were significantly enriched in catalytic activity and binding.

KEGG Pathway Classification of DEGs
Since there is no reference genomic sequence for fat greenling, the function and classification of the selected genes are obtained by comparison with sequence homology analysis of related gene sequences of the species in Table 5. Function and pathway allocation of DEGs were performed after ammonia exposure 48 and 96 h, and all DEGs were analyzed against the KEGG database. The top 20 pathways with the most abundant DEGs involved in with ammonia exposure are shown in Figure 8, and pathway-enrichment analysis showed that DEGs were mainly enriched in carbon metabolism and biosynthesis of amino acids. After ammonia exposure for 48 and 96 h, DEGs enriched in carbon metabolism were dominant. DEGs induced FIGURE 4 | GO classification of assembled unigenes. X-axis: GO classification, including biological process, cellular component, and molecular function, Y -axis: The left ordinate is the percentage of genes; the right ordinate is the numbers of genes.  by exposure mainly were reflected in environmental information processing that play an important role in immune function, oxidative stress, and apoptosis. Seven subclasses related to immune function, oxidative stress, and apoptosis, and further elaborate on the basis of KEGG pathway analysis, including the MAPK-signaling pathway, FOXO-signaling pathway, phosphatidylinositol-signaling pathway, cytokine-cytikine receptor interaction, Hippo-signaling pathway, and neuroactive ligand-receptor interactions ( Table 5). A total of 16 DEGs were induced, of which three genes were up-regulated and 13 genes were down-regulated. These DEGs are associated with oxidative stress and apoptosis.

Validation of Expression Genes by qRT-PCR
To verify RNA-seq results, we randomly selected four upregulated and six down-regulated DEGs with different fold changes for qRT-PCR analysis. As shown in Figure 9, in both experiments, the expression of these 10 genes showed the same pattern of up-regulation and down-regulation. The qRT-PCR analysis results are thus consistent with data from the RNA-seq.

DISCUSSION
Fat greenling is one of the most common commercial aquatic animals. Studies on fat greenling mainly have focused on pollution biology and physiology; few have focused on transcriptome analysis, with the exception of spleen analysis following V. harveyi infection (Diao et al., 2019). We carried out high-throughput sequencing nine samples of fat greenling under ammonia exposure. A total of 78.48 Gb of clean data were obtained, the clean data of each sample reached 6.02 Gb, and the Q30 was above 89.97%, indicating that the sequencing effect was good and provided a credible basis for subsequent assembly. A total of 87,675 unigenes were obtained after assembly, of which 16,756 unigenes were over 1 kb in length, 28,062 unigenes were successfully annotated by blast alignment databases. The experimental sequencing data coverage was high, and the library construction was successful. These results contribute to our understanding of fat greenling genes, providing data for the future study of the species.
Ammonia is a highly toxic environmental pollutant (Ye et al., 2018). Previous reports have suggested that high concentrations of ammonia in water accumulate in different tissues of aquatic organisms, triggering oxidative stress, immunosuppression, and other physiological problems (Cheng et al., 2015;Zhang et al., 2015;Li et al., 2016). GO-enrichment analysis showed that DEGs in fat greenling liver under ammonia exposure were significantly enriched in terms of cellular process, metabolic process, and single-organism process in BF, which were correlated with oxidative stress and immunological function. KEGG analysis showed that the DEGs were mainly involved in environmental information processing, including the MAPK-signaling pathway, FOXO-signaling pathway, phosphatidylinositol-signaling pathway, cytokine-cytikine receptor interaction, Hippo-signaling pathway-fly, and neuroactive ligand-receptor interactions. These pathways were mainly associated with oxidative stress and immunological functioning. Several important KEGG pathways to be discussed below are likely to be involved in the response of fat greenling under ammonia exposure.

MAPK Signaling Pathway
Eukaryotic cells regulate gene expression through a cascade of MAPK (mitogen-activated protein kinase) signals under various stressful conditions, such as oxidative stress inflammation (Obata et al., 2000;Kyriakis and Avruch, 2001;Lee and Choi, 2003;Marshall et al., 2005;Feidantsis et al., 2012). In this study, some DEGs in the MAPK pathway were associated with cell damage, disease, and immunity, including NF-kappa-B, fibroblast growth factor (FGF-13), mitogen-activated protein kinase kinase kinase 8 (MAP3K8), interleukin-1 receptor (IL-1RI), proto-oncogene and dual-specificity protein phosphatase 2 (Dusp2) ( Table 5). NF-kappa-B is a pleiotropic transcription factor existing in almost all cell types and is the end point of a series of signal transduction events triggered by a series of stimuli related to biological processes such as inflammation, immunity, differentiation, cell growth, tumorigenesis, and apoptosis (Lee et al., 2013). Studies have shown that medaka defensins are important cationic  antimicrobial peptides, which have been shown to have a relationship with NF-kappa-B and play an important role in resistance to pathogens (Zhao, 2008). Down-regulation of NF-KB is triggered by autoimmune diseases, chronic inflammation, and cancers. MAP3K8 (TPL-2) is associated with the TLR-signaling pathway; TLRs target a common signaling pathway that causes activation of NF-KB, MAPKs, ERK, P38, and JNK. When cells receive signal stimulation, p105 is phosphorylated (Salmeron et al., 2001), causing its ubiquitination to release tpl2; the released tpl2 is activated by phosphorylation of the site, thereby activating MEK1/2 and its downstream ERK1/2, realizing the entire MAPK signal (Lang et al., 2003).
Interleukin type I receptor (IL-1RI) intracellular regions are also present in the Drosophila Toll protein, a discovery that highlights the role members of the il-1 family have in innate responses (Gay and Keith, 1991). IL-1β is a major inflammatory cytokine that also plays a role in inflammation associated with many chronic diseases (Wettlaufer et al., 2016). The expression of IL-1RI in kidney, spleen, liver, and sputum of Salmo salar was up-regulated after LPS and TNF-α stimulation, indicating that IL-1RI plays a major role in inflammatory response by regulating IL-1β (Subramaniam et al., 2002).
Under normal conditions, intracellular c-fos is highly conserved and has only low levels of expression which are difficult to detect (Blume et al., 1999). When the body is stimulated by external factors, intracellular c-fos can be induced rapidly (Usuki et al., 2000;Dong et al., 2005). A variety of extracellular signals can activate downstream transcription factors through second messengers, induce Fos protein biological expression, act on target genes, and then respond to external stimuli (Akins et al., 1996). At present, a variety of fish c-fos genes have been cloned, including in Oncorhynchus mykiss (Matsuoka et al., 1998), Rivulus marmoratus (Li et al., 2004), Dicentrarchus labrax (Rimoldi et al., 2009), Fugu rubripes (Trower et al., 1996), and Oryzias latipes (Okuyama et al., 2011). We found that ammonia exposure led to the expression of the MAPK-signaling pathway in the liver of  fat greenling, and the response is used to mobilize the upregulated and down-regulated expression of related genes to cope with stress.

FOXO-Signaling Pathway and Phosphatidylinositol-Signaling Pathway
Studies have shown expression of related factors in the FOXO-signaling pathway, which can increase expression of related antioxidant defense system enzymes (Cox and Lane, 1995). Down-regulated expression of phosphoinositide 3-kinase (PI-3K) genes occurs in both the FOXO-signaling pathway and phosphatidylinositol-signaling pathway. Phosphoinositide 3-kinase is an enzyme with serine/threonine protein kinase activity, consisting of two chains, P85 and P110, and is involved in immune cell activation signal transduction and activation. Reports have shown that in T cells, a variety of cytokines, antigen receptors, and environmental factors stimulate activation of PI-3K. Akt is activated accordingly, as well as other downstream transcription factors, in response to viruses and other environmental factors. As such, they play a key role in maintaining a positive immune response of T cells, and promote development of T cells during the immune response (Sasaki et al., 2002). We showed expression of G1/s-specific cyclin-d1 and the PI-3K gene were down-regulated, and we speculate that ammonia had a negative regulatory effect on the PI-3K signal pathway. We observed that, compared with the control group, fish showed sluggish activity, and some fish had ulcerations on their skin, indicating that stress resistance of fish may have been decreased.

Cytokine-Cytikine Receptor Interaction and Hippo-Signaling Pathway-Fly
Chemokines are low-molecular (mostly 8-10KD) proteins that are capable of attracting leukocytes to the site of infection, and they play an important role in the inflammatory response. C-X-C chemokine receptor type 3 (CXCR3) plays an important role in regulating the transport and function maintenance of T cells (Groom, 2011). There are many studies that have focused on induction of chemokines and chemokine receptors by pathogens and environmental stimuli. For example, acute toxicity of sodium nitrite, along with bacterial, fungal, and viral infections, can induce the expression of CXCR3a in Channa striatus (Bhatt et al., 2014). LPS also can induce the expression of CXC chemokines in carp (Savan et al., 2003).
We found CXCR3 involved in cytokine and cytokine receptor pathways was up-regulated, as well as were genes in other immune-related pathways, including interleukin-1 receptor type 2-like, fibroblast growth factor 13-like isoform X2, protooncogene c-fos like, and dual-specificity protein phosphatase 2-like. These are helpful for screening immune-related genes and analyzing their expression.
At present, there is no report thorough published examination on the HIPPO signaling pathway in fishes. YAP1 also is associated with apoptosis. It is speculated that the down-regulation of YAP1 gene under ammonia-N stress is related to the occurrence of various diseases.

Neuroactive Ligand-Receptor Interaction and Hippo Signaling Pathway-Multiple Species
Prostaglandin E2 (PGE2) is an important arachidonic acid metabolite with strong immune activity, which can regulate the development, function, and survival of immune cells (Tilley et al., 2001). E-prostanoid receptors (EPs) are divided into four subtypes: EP1, EP2, EP3, and EP4. PGE2 ACTS on T lymphocytes and dendritic cells through EP4 receptors to promote differentiation of Th1 cells and the amplification of Th17 cells. Recent studies have shown that PGE2-EP4 signaling plays an important role in some autoimmune diseases (Libioulle et al., 2007;Yao et al., 2009;Chen et al., 2010;Ngoc et al., 2011).
In summary, our results show that fat greenling experience ammonia-driven physiological stress. A number of DEGs were significantly enriched, and related to pathways of immune or oxidative stress. Parallel results from GO functional analysis and KEGG-signaling pathway classification analysis, suggest we have identified various genes related to immune and oxidative stress. Ammonia exposure clearly seems to affect liver function, and our data provide a foundation for future research on this ecologically and economically important fish species.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI Sequence Read Archive Database under accession number SAMN11155636.

ETHICS STATEMENT
The animal study was reviewed and approved by the Marine Biology Institute of Shandong Province.

AUTHOR CONTRIBUTIONS
FH and WG: substantial contributions to the conception and design of the work and provide approval for publication of the content. LL, FG, YJ, XuW, XiW, LP, and DL: acquisition of data for the work. LL: analysis and interpretation of data for the work. LL and FH: drafting the work and revising it critically for important intellectual content. All authors agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy and integrity of any part of the work are appropriately investigated and resolved.

FUNDING
This study was supported by the Key Technology Research and Development Program of Shandong (2019GHY112071 and 2019GHY112062) and Major agricultural application technology innovation projects in Shandong Province (SD2019YY007).