Time-series transcriptomic analysis of the kelp grouper Epinephelus moara in response to low salinity stress

ABSTRACT The Kelp grouper Epinephelus moara is one of the most widely consumed and economically important marine fish in China. The species can tolerate a wide range of salinity, but genomic resources are not available, and the molecular mechanisms underlying adaptation to salinity at the transcriptomic level remain largely unclear. In this study, the transcriptomic responses of the liver of E. moara under low salinity were investigated using the Illumina digital gene expression system. After de novo assembly, 499,356 transcripts were generated and contributed 445,068 unigenes. A total of 14, 19, 33 and 3101 genes were differentially expressed following exposure to low salinity stress for 2, 6, 24 and 48 h, respectively. Only two genes were differentially expressed in all groups. Four genes related to metabolism and ambient salinity adaption were randomly selected to validate the differentially expressed genes (DEGs) by real-time PCR. Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analysis were used to analyse the functional significance of DEGs, including those responding to salinity through diverse biological processes, cellular components, molecular functions, and pathways associated with metabolic and osmotic responses. This work provides new insight into the response to salinity challenges in E. moara, and the findings expand our knowledge of the molecular basis of metabolic regulation mechanisms in this species. Additionally, the transcriptional data provide a valuable resource for future molecular and genetic studies on E. moara.


Introduction
The Kelp grouper Epinephelus moara is a large marine fish that inhabits the bottom of warm oceans in the western North Pacific, particularly those surrounding Vietnam, southern and eastern China, South Korea and Japan (Guo et al. 2009). The meat of this species is renowned for its delicious, fresh taste, and it is nutrient-rich, making it a favourite of consumers and bestowing high economic value. Due to the limited abundance of E. moara populations and overfishing, the species is gradually becoming scarce and the price keeps rising. In China, many researchers have carried out artificial breeding and seed cultivation experiments on E. moara, and found that it grows fast and possesses high adaptability (Tian et al. 2017a). The fish is an important marine aquaculture species in south-eastern coastal regions of China.
Salinity is an environmental factor affecting osmotic pressure regulation in fish that directly influences metabolism and immune responses, restricting the distribution of fish in natural ecosystems. In mariculture, salinity has a significant effect on the development of fish embryos, larvae and juveniles (Sørensen et al. 2016;Kujawa et al. 2017). It affects the ability of grouper to regulate osmotic pressure, the activity of metabolic enzymes, and the abundance of immunerelated factors, influencing the ability to regulate osmotic pressure, digestion and absorption of nutrients, resistance to disease, and ultimately growth and survival (Tsui et al. 2012;Cheng et al. 2013;García et al. 2013;Chapman et al. 2014;Sutthinon et al. 2015;Chen et al. 2016;Arrokhman et al. 2017). Therefore, studying the effects of salinity on the growth and development of grouper is important for artificial breeding programmes. Due to the widely fluctuating salinity in coastal areas of China, especially during the breeding season when rainfall can be high and salinity low in estuaries, the growth and survival of aquatic animals can be restricted. Therefore, studying the impact of low salinity stress on E. moara is important.
Salinity is an important factor in the successful recruitment of larval and juvenile groupers (Gracia-López et al. 2004;Cheng et al. 2013). Various enzymes, transporters and signal transmission pathways participate in the processes of salinity adaption and osmoregulation to maintain an internal osmotic and ionic homeostasis when live in a salinity changing environment (Tsui et al. 2012;Sutthinon et al. 2015;Chen et al. 2016). It is believed that the drastic change of liver metabolism can be triggered by changing the salinity . Analysis and identification of candidate genes involved in salinity change is the first step to elucidate the molecular mechanism and understand factors underlying this core physiological process. Transcriptomics can be used to investigate the transcriptional regulation of all genes in cells and tissues. The development and application of transcriptomics approaches can not only reveal details of the cellular components, but also help to elucidate the functional genes and components in the genome, which are of great importance for studying the growth and immune responses of organisms (Meng et al. 2015;Cao et al. 2017;Zhang et al. 2017). In recent years, transcriptomics has been applied by aquaculture researchers to study the effects of salinity stress on aquatic animals (Thanh et al. 2014). In the present study, next-generation high-throughput sequencing was used to analyse key genes and signal transduction pathways related to salinity stress in E. moara. The findings deepen our understanding and will be of benefit to grouper breeding programmes and controlling diseases during aquaculture of this important food species.

Materials and methods
Experimental fish E. moara were obtained from the Yellow Sea Fisheries Research Institute. This study was conducted at the Shanghai Fisheries Research Institute. Fish with a similar body length (10.43 ± 1.08 cm) and weight (16.73 ± 4.47 g) were acclimatised in aerated normal seawater with a salinity of 27% for 2 weeks before reduced salinity exposure, during which 25% of the tank water was renewed daily. After acclimatisation, 60 healthy fish were randomly divided into three tanks and exposed to a low salinity of 15 ppt. Three individual fish from each tank were randomly sampled after exposure to the reduced salinity conditions for 0, 2, 6, 24, and 48 h. Liver samples were collected and immediately flash-frozen in liquid nitrogen and stored at -80°C until they were processed for RNA extraction.

RNA extraction and library preparation for Illumina sequencing
Total RNA from livers was extracted using Trizol reagent (Invitrogen, USA) according to the manufacturer's instructions. RNA quantity, purity and integrity were determined by 1% agarose gels and a Nanodrop 1000 spectrophotometer (NanoDrop, USA). High-quality RNA was selected for high-throughput sequencing, and samples from the three fish for each time point were combined. mRNA was purified from total RNA using oligo-dT-attached magnetic beads and fragmented using divalent cations under elevated temperature. First-strand cDNA was synthesised using random primers and reverse transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. A pairedend library was constructed using the Genomic Sample Prep Kit and sequenced using the Illumina Mi-Seq platform (Yin et al. 2016).

Assembly and annotation of the transcriptome
Raw sequencing reads with an average length of 250 bp were generated and filtered through the Illumina sequencing platform. All dirty raw reads with adaptors, unknown bases at >10%, and low-quality sequences were trimmed to generate a set of clean reads. Short sequences were also removed using a custom Perl Program. The resulting high-quality sequences stemming from quality control analysis were de novo assembled into contigs and transcripts with Trinity software. To reduce data redundancy and evaluate the efficacy of the assembly, transcripts with a minimum length of 200 bp were assembled and clustered using TGICL under default parameters. The longest transcript among the multitude of assembly transcript isomers was defined as a unigene. All unigenes were used as queries for searching nr, nt and Swiss-Prot databases, and functionally annotated by Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis (Song et al. 2016).

Identification of differentially expressed genes (DEGs)
The number of reads aligned to each unigene in the alignment file and the read counts were normalised as reads per kilobase per million (RPKM) values. The false discovery rate (FRD) was used to determine the threshold of the p-value in multiple tests to judge the significance of gene expression differences. DEGs between pairwise comparisons (0 h vs. 2 h, 0 h vs. 6 h, 0 h vs. 24 h, and 0 h vs. 48 h) were identified using the Bioconductor tool edgeR (p-value <0.05, FDR <0.05). Functional enrichment analysis of DEGs was performed by mapping each DEG to GO and KEGG databases. The FDR was used to correct the p-value and, pathways with FDR values <0.05 were considered significantly enriched (Thanh et al. 2014;Nie et al. 2016).

Quantitative real-time PCR (qRT-PCR) analysis
To validate our Illumina sequencing data, four genes were randomly selected for qRT-PCR analysis using a SYBR Premix Ex Taq kit (Invitrogen) according to the manufacturer's instructions. The same pooled RNA samples used for generating RNA-seq data were reverse-transcribed into cDNA using a PrimeScript 1 st strand cDNA synthesis kit (Takara, Japan). Specific primers used for qRT-PCR were designed using Primer 5 software and are listed in Table 1. The β-actin gene was chosen as an endogenous control. Relative gene expression was determined using the 2 −△△Ct method Yin et al. 2016).

Transcriptome sequencing and de novo assembly
A total of 63,545,722, 37,741,048, 49,822,946, 64,207,362 and 126,858,410 raw reads were generated from the 0, 2, 6, 24 and 48 h groups, respectively. After removing reads containing adapter, reads containing poly-N sequences, and low-quality reads, 61,757,036, 36,311,164, 47,969,516, 62,429,208 and 124,017,584 clean reads with an average length of 140, 139, 142, 141 and 138 bp were obtained for subsequent analysis. The percentages of aligned reads were 57.3%, 55.9%, 66.9%, 60.4% and 51.7%, respectively. Further assembly analysis showed that these high-quality reads generated 499,356 transcripts and contributed to 445,068 unigenes. The quality of these reads is shown in Tables 2 and 3.

Identification and analysis of DEGs
Overall, 3294 genes from the four time point groups were significantly differentially expressed in response to salinity acclimation. Of these, 3167 were up-regulated significantly (p-value <0.05), comprising 14 (0 h vs. 2 h), 19 (0 h vs. 6 h), 33 (0 h vs. 24 h) and 3101 (0 h vs. 48 h) genes in the four respective groups. The number of genes down-regulated (p-value <0.05) was 127, comprising 21, 29, 37 and 40 genes in the 2, 6, 24 and 48 h groups (Table 4). However, only two DEGs were present in all pairwise comparisons (Figure 1(a)). Many more DEGs (3167) were up-regulated between 0 and 48 h (Table 4), and a wide variety of genes were up-regulated following salinity stress (Figure 1(b)). TRAF-type zinc finger domain-containing protein 1-like  Table 4). The cellular component category was highly represented by MHC class II protein, complex MHC protein and endosome sub-categories. The biological process category was highly represented by peptide antigen binding, phosphate ion binding, and antigen binding sub-categories, while antigen processing and presentation of exogenous peptide antigen via MHC class II, peptide antigen via MHC class II, and peptide or polysaccharide antigen via MHC class II were the dominant biological process subcategories (Figure 2(a)).
In the 0 and 6 h pairwise comparison, 265 GO terms were significant (p-value <0.05), including 38 molecular function, 42 cellular component, and 185 biological process sub-categories (Table 4). In the 0 and 24 h pairwise comparison, 308 GO terms were identified (pvalue <0.05), including 184 biological process, 57 cellular component, and 67 molecular function sub-categories (Table 4). Furthermore, in the 0 and 48 h pairwise comparison, the number of GO terms (p-value <0.05) was increased significantly (2155), including 1427 biological process, 188 cellular component, and 540 molecular function sub-categories (Table 4). In general, GO terms identified in the 0 and 48 h pairwise comparison were much more significant than from other pairwise comparisons. The top 10 most enriched GO terms in the three main GO categories (biological process, cellular component and molecular function) among all pairwise comparisons are shown in Figure 2Figure 2(a, c, e and g).

KEGG analysis
KEGG pathway analysis was performed to investigate pathways that were significantly altered following exposure to low salinity. The most enriched KEGG pathways are shown in Figure 2

Validation of RNA-seq data by real-time PCR
To validate the sequencing data, four unigenes, namely TRAF-type zinc finger domain-containing protein 1-like, parvalbumin beta, SH3 domain-binding glutamic acidrich-like protein 3-like, and cold shock domain-containing protein E1, were randomly selected for real-time PCR experiments using the same pooled RNA samples employed for generating the RNA-seq data. The expression levels of these four genes measured by RNA-seq and real-time PCR are presented in Figure 3. The expression levels of all four genes measured by real-time PCR were consistent with those of RNA-seq, confirming the validity of the RNA-seq results.

Discussion
Salinity is an important factor in the marine environment that can have a direct impact on the growth and development, reproduction, and immunity of fish and other Changes in salinity can force fish to adjust the dynamic balance of the osmotic pressure through a series of physiological changes, which will lead to corresponding changes physiological and biochemical processes such as digestion and absorption of food, nutrient metabolism, and immunoregulation (Tran-Ngoc et al. 2017;Tian et al. 2017b). Salinity changes can directly affect the activity of enzymes regulating osmotic pressure and gene expression. Therefore, studies on the effects of salinity on fish have mainly focussed on osmoregulatory organs (gill, intestine and kidney) (Havird et al. 2016;Zhang et al. 2016). Although the liver is not an osmoregulatory organ, it is the major source of energy for osmoregulatory organs. In fish, the liver is important for the metabolism of sugars, fats, proteins, vitamins, hormones, and other components, as well as detoxification, secretion of bile, immune defences, and coagulation, hence the intense focus on this organ in studies on fish metabolism (Cuesta et al. 2005;Singh et al. 2017;Wu et al. 2017). Although energy metabolism in the liver is likely enhanced during salinity adaptation, little is known about the reorganisation of metabolism in this organ during this process in fish. Previous studies on osmoregulation and liver metabolism focused on physiological processes related to osmotic pressure, and energy production and expenditure, but gene expression profiling is not widely reported. Therefore, the present study focused on expression profiling of the liver in grouper exposed to low salinity stress. Transcriptomic analysis can probe gene expression, genetic structures and functions at a genome-wide scale, which is likely to be useful for elucidating and analysing the mechanisms underlying salinity stress adaptation in E. moara.
In recent years, transcriptomic analysis of fish in response to salinity stress has mainly focused on genes and enzymes related to osmotic organs and osmotic pressure regulation (Havird et al. 2016), and studies on liver metabolism have not yet been reported. The effects of low salinity stress on gene expression in the liver of E. moara were investigated for the first time in this study using the Illumina digital gene expression system. The results showed that the number of DEGs increased with increasing duration of low salinity stress, suggesting this species constantly adjusts in order to adapt to low-salinity environments. Within the first 24 h of low-salinity stress, although the number of DEGs increased continuously, the total number was only 35, 48 and 70 at 2, 6 and 24 h, respectively. By contrast, after low-salinity stress for 48 h, 3141 DEGs were Figure 3. Validation of RNA-seq data by real-time PCR. The expression of TRAF-type zinc finger domain-containing protein 1-like, parvalbumin beta, SH3 domain-binding glutamic acid-rich-like protein 3-like, and cold shock domain-containing protein E1 was measured by RNA-seq (white columns) and real-time PCR (black columns) and plotted as -log2(fold change).
identified, which indicates that sustained low-salinity stress has a significant effect on metabolism in fish. In contrast to the results of this study, Jian et al. (2017) found that the number of DEGs initially increased then decreased with prolonged salinity stress (30, 60 and 90 min). Other studies (Gibbons et al. 2017;Jian et al. 2017) found that whether comparing between salinity levels or between salinity stress duration, only a relatively small number of DEGs were identified, consistent with the results of the present study, in which only 2 DEGs were present in all time point groups.
Real-time PCR is widely used to indicate the reliability and accuracy of RNA-seq data, and while both techniques provide information on gene expression levels, analysis of the functions of DEGs and their associated pathways can be fruitful. GO can be used to define and characterise the functions of DEGs in terms of biological process and cellular components based on the relevant molecular functions of the gene products (Meng et al. 2015;Nie et al. 2016). GO analysis involves annotating unknown genes with known functions from similar genes in other species based on molecular function, biological process and cellular component categories. In the present study, the number of GO terms belonging to the biological process was much larger than those belonging to the molecular function and cellular component categories. With increasing low-salinity stress duration, the number of GO terms in all three categories was increased, and the significance (-log2[p-value]) was also increased. Among the terms, those related to plastid and chloroplast cellular component sub-categories were the most abundant. These results suggest the metabolism in the liver of E. moara was significantly affected by the duration of low-salinity stress, leading to changes in liver cell components.
In addition to GO analysis, KEGG pathway analysis was performed to elucidate the probable functional status of assembled transcripts. KEGG analysis can provide a more comprehensive and complex analysis of cellular metabolic pathways and processes such as signal transduction, enzymatic reactions, and membrane transport, potentially providing greater biological information . KEGG pathway analysis has been widely used to study the effects of salinity stress on fish, but most studies have focused on osmotic pressure regulation (Thanh et al. 2014;Zhang et al. 2016), while liver metabolism has been largely ignored. Our results showed that within the first 24 h of low-salinity stress, KEG pathways were not significantly altered (-log2[pvalue]). However, after 48 h of low-salinity stress, KEGG pathways were altered significantly, and metabolic pathways were the most different, indicating damage to liver tissue following prolonged exposure to low-salinity conditions.
In conclusion, the results of this study provide a first glimpse into the changes in global gene expression in the liver of E. moara that occur in response to low salt stress. A total of 35, 48, 70 and 3141 DEGs were identified after 2, 6, 24 and 48 h of exposure to low-salinity conditions. GO and KEGG enrichment analyses revealed that the expression of genes associated with metabolism were influenced most significantly by low-salinity stress. These pathways and gene expression patterns expand our understanding of the molecular mechanisms underpinning the metabolic changes that occur in response to salinity stress in this important aquaculture species.