RNA Sequencing Analysis Reveals Divergent Adaptive Response to Hypo- and Hyper-Salinity in Greater Amberjack (Seriola dumerili) Juveniles

Simple Summary The gill tanscriptomes of greater amberjack (Seriola dumerili) reared under different salinity stress were analyzed. The regulatory networks of salinity-related pathways were explored through Kyoto Encyclopedia of Gene and Genome (KEGG) pathway enrichment and bioinformatics analyses. This will be of great value in understanding the molecular basis of salinity adaptation in greater amberjack. Abstract Salinity significantly affects physiological and metabolic activities, breeding, development, survival, and growth of marine fish. The greater amberjack (Seriola dumerili) is a fast-growing species that has immensely contributed to global aquaculture diversification. However, the tolerance, adaptation, and molecular responses of greater amberjack to salinity are unclear. This study reared greater amberjack juveniles under different salinity stresses (40, 30, 20, and 10 ppt) for 30 days to assess their tolerance, adaptation, and molecular responses to salinity. RNA sequencing analysis of gill tissue was used to identify genes and biological processes involved in greater amberjack response to salinity stress at 40, 30, and 20 ppt. Eighteen differentially expressed genes (DEGs) (nine upregulated and nine downregulated) were identified in the 40 vs. 30 ppt group. Moreover, 417 DEGs (205 up-regulated and 212 down-regulated) were identified in the 20 vs. 30 ppt group. qPCR and transcriptomic analysis indicated that salinity stress affected the expression of genes involved in steroid biosynthesis (ebp, sqle, lss, dhcr7, dhcr24, and cyp51a1), lipid metabolism (msmo1, nsdhl, ogdh, and edar), ion transporters (slc25a48, slc37a4, slc44a4, and apq4), and immune response (wnt4 and tlr5). Furthermore, KEGG pathway enrichment analysis showed that the DEGs were enriched in steroid biosynthesis, lipids metabolism, cytokine–cytokine receptor interaction, tryptophan metabolism, and insulin signaling pathway. Therefore, this study provides insights into the molecular mechanisms of marine fish adaptation to salinity.


Introduction
Marine environmental factors, such as salinity, low O 2 concentration, temperature, and pH value, influence the physiological and biological status of marine animals [1]. For instance, environmental stresses activate the sympathetic nervous system [2], the release of adrenaline and noradrenaline [3], and the hypothalamic-pituitary-interrenal axis in fish [4],

Experimental Fish, Salinity Development, and Tissue Collection
A total of 80 greater amberjack juveniles (body length, 8.33 ± 0.45 cm and body weight, 6.38 ± 1.33 g) were used in this study. Before experiments, the greater amberjack juveniles were reared in tanks at 22 ± 1.0 • C in Donghai Island (Zhanjiang, Guangdong, China). They were randomly divided into four cylindrical 1000 L tanks (20 individuals per tank) at various salinities: 40, 30, 20, and 10 ppt (parts per thousand) groups. The salinity levels in the experiment were selected based on the previous salinity adaption study of greater amberjack larvals [47] and our unpublished data of juveniles. The fish in the control group were reared in natural seawater with a salinity of 30 ppt. The 40, 20, and 10 ppt groups were regulated using a commercial seawater salty crystal and aerated tap water. The fish were fed on commercial float bait twice a day at 9:00 and 19:00 for 30 days. Six fish were randomly selected from each group on 0, 15, and 30 days in our study, which was based on studies of other fishes, such as Asian seabass (Lates calcarifer, Bloch, 1790) [48], catfish (Lophiosilurus alexandri) [24], and cobia (Rachycentron canadum), under different salinity stress levels [31]. Unfortunately, the fish in the 10 ppt group were all dead within 10 days. The fish were then anesthetized using 100 mg/L tricaine methane sulfonate (MS 222; Sigma-Aldrich, St. Louis, MO, USA) and dissected. The gill tissues were immediately collected in centrifuge tubes containing 1 mL RNA stabilization reagent overnight, then stored at −80 • C for RNA extraction, sequencing, and gene expression analysis. The RNA of gill samples after 30 days were used for transcriptomic analysis, all gill samples after 15 and 30 days were used for qPCR verification.

Total RNA Extraction, Library Construction, and Illumina Sequencing
Trizol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from the gill tissues, following the manufacturer's instructions. As previously described, the cleavage of gill tissue samples, RNA extraction, RNA purity, degradation, and contamination examinations were performed [49]. An Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) was used to detect RNA integrity. Total RNA with an RNA integrity number (RIN) score >7 was used for sequencing.
TrueSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA) was used to obtain complementary DNA (cDNA) libraries from gill tissue, following the manufacturer's instructions. Then, 3 µL of USER Enzyme (NEB, New England Biolabs, Palo Alto, CA, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min, followed by 5 min at 95 • C before PCR. PCR was performed using Phusion High-Fidelity DNA polymerase, universal PCR primers, and Index (X) Primer. The AMPure XP system was used to purify the PCR products. An Agilent Bioanalyzer 2100 system was used to assess the library quality. A cBot Cluster Generation System was used to cluster the index-coded samples via the TruSeq PE Cluster Kit v4-cBot-HS (Illumia, San Diego, CA, USA), following the manufacturer's instructions. The library preparations were then sequenced on a HiSeq X-ten platform to generate the paired-end reads (150 bp).

Transcriptome Assembly and Functional Gene Annotation
The assembled Seriola dumerili genome (deposited in the DNA Data Bank of Japan (DDBJ) under accession numbers BDQW01000001-BDQW01034655 (Biosample ID: SAMD00083043_sdu_WGS.acclist.zip)) was used as a reference database for mapping reads [50]. The Illumina high-throughput sequencing platform was used to sequence the cDNA library, generating raw reads/data. The raw reads/data were filtered, then the adapter sequence and ploy-N (unable to determine base information) and low-quality (reads with <50% bases of quality value) reads were removed to obtain high-quality clean reads/data. Raw data and clean data were saved in FASTQ format. The Q20, Q30, GCcontent, and sequence duplication levels of the clean data were measured. All the subsequent analyses were based on high-quality clean data. BLASTx

Analysis of Differential Expressed Genes (DEGs) Analysis
The greater amberjack gene expression levels were detected using the fragments per kilobase per million (FPKM) method. The DEGs in 20 vs. 30 ppt and 30 vs. 40 ppt groups were identified using the DESeq2 R package (version 1.16.1). DESeq2 was used to determine differential expression in digital gene expression data based on the negative binomial distribution model. Genes with a fold change ≥2 and a false discovery rate (FDR) < 0.05 were considered DEGs. A KEGG pathway analysis was conducted, and DEGs with p < 0.05 were considered statistically significant [31]. In Donghai island, 30 ppt was the natural salinity of seawater. The DEGs in both 20 vs. 30 ppt and 30 vs. 40 ppt were identified and screened in order to clarify the effects of hypo-salinity and hyper-salinity seawater environments on gills of greater amberjack compared with the natural salinity.

Quantitative Real-Time PCR (QPCR) Validation
The expression patterns of DEGs in the RNA sequencing analysis were validated using qPCR. Sample collection, RNA extraction, and reverse transcription were performed as previously described [51]. RNA samples were obtained from the 20, 30, and 40 ppt groups. Each group had three replicates. The DEGs primers were designed using Primer5 software based on the assembled transcripts. A light CyclerTM96 (Roche, Indianapolis, IN, USA) was used for qPCR analysis, following the protocol of SYBR Green Real-Time PCR Master Mix (Takara, Tokyo, Japan). β-actin was used as a reference gene to normalize the expression levels [52][53][54]. The relative abundance of DEG mRNA transcripts was evaluated using the 2 − Ct method. The primer sequences for qPCR are shown in Table S1.

Statistical Analysis
The relative mRNA expression levels and FPKM values are expressed as mean ± standard error (SE). The one-way ANOVA with Tukey's post hoc test was used to evaluate the significant differences between 20 vs. 30 ppt and 30 vs. 40 ppt groups. The significance level was set at α = 0.05. The Statistical Package for the Social Sciences (SPSS) 16.0 (SPSS, Chicago, IL, USA) was used for all statistical analyses.

Illumina Sequencing
The gill transcriptome can provide a valuable RNA resource for future analysis of greater amberjack adaptability to salinity and artificial culture. The HiSeq X Ten platform was used for RNA sequencing of gill samples at 20, 30, and 40 ppt. A total of 94.79, 93.66, and 78.83 million clean reads were obtained at 20 (G20), 30 (G30), and 40 ppt (G40), respectively, after quality control. The Q30 values and GC content of the clean reads were more than 95% and 49%, respectively. The high-quality reads were used for further analysis. The reads' Q20 and Q30 values, GC content, and transcript numbers for each cDNA library are shown in Table 1.  Table 2). The annotated genes provided the basis for further analysis of the specific molecular processes in greater amberjack. The best BLAST results of reads were enriched for closely related fish species, including Seriola lalandi (6.24%), Lates calcarifer (2.09%), Larimichthys crocea (1.49%), Stegastes partitus (0.63%), Oreochromis niloticus (0.55%), and other species (4.98%) ( Figure 1).

Identification and Analysis of Differentially Expressed Genes (DEGs)
A total of 417 (205 up-regulated and 212 down-regulated) and 18 (nine up-regulated and nine down-regulated) DEGs were identified in the G30 vs. G20 and G30 vs. G40 groups, respectively, using DESeq2 software, FDR-adjusted p-value < 0.05 and |Log2(fold change)| ≥ 1 ( Figure 2). Heat maps of the clustered DEGs under hypo-and hyper-salinity stresses are shown in Figure 3. The top 20 DEGs in the G30 vs. G20 group and the top 18 DEGs in the G30 vs. G40 group are shown in Table 3.

Identification and Analysis of Differentially Expressed Genes (DEGs)
A total of 417 (205 up-regulated and 212 down-regulated) and 18 (nine up-regulate and nine down-regulated) DEGs were identified in the G30 vs. G20 and G30 vs. G4 groups, respectively, using DESeq2 software, FDR-adjusted p-value < 0.05 and |Log2(fol change)| ≥1 ( Figure 2). Heat maps of the clustered DEGs under hypo-and hyper-salinit stresses are shown in Figure 3. The top 20 DEGs in the G30 vs. G20 group and the top 1 DEGs in the G30 vs. G40 group are shown in Table 3.

Identification and Analysis of Differentially Expressed Genes (DEGs)
A total of 417 (205 up-regulated and 212 down-regulated) and 18 (nine up-regulated and nine down-regulated) DEGs were identified in the G30 vs. G20 and G30 vs. G40 groups, respectively, using DESeq2 software, FDR-adjusted p-value < 0.05 and |Log2(fold change)| ≥1 ( Figure 2). Heat maps of the clustered DEGs under hypo-and hyper-salinity stresses are shown in Figure 3. The top 20 DEGs in the G30 vs. G20 group and the top 18 DEGs in the G30 vs. G40 group are shown in Table 3.

DEG Annotation and Pathway Analysis
Some enriched KEGG pathways related to Seriola dumerili metabolism and molecular signaling pathway were identified, including steroid biosynthesis, cytokine-cytokine receptor interaction, porphyrin and chlorophyll metabolism, cytosolic DNA-sensing pathway, the intestinal immune network for IgA production, and tryptophan metabolism (Fig-Figure 3. Expression heatmap of the DEGs in G30 vs. G20 group (a) and G30 vs. G40 group (b). Green and red squares indicate down-and up-regulation, respectively. The brighter colors indicate the significant fold changes.

DEG Annotation and Pathway Analysis
Some enriched KEGG pathways related to Seriola dumerili metabolism and molecular signaling pathway were identified, including steroid biosynthesis, cytokine-cytokine receptor interaction, porphyrin and chlorophyll metabolism, cytosolic DNA-sensing pathway, the intestinal immune network for IgA production, and tryptophan metabolism (Figures 4 and 5). The top 10 enriched KEGG pathways of the DEGs under hypo-and hyper-salinity stresses are shown in Table 4.

Validation of RNA Sequencing Data by QPCR
Considering the top 20 most important salinity stress responses of DEGs identified in the transcriptomic analysis, which had been reported in the previous studies, nine genes were selected for RNA sequencing data validation using qPCR ( Figure 6). The mRNA expression of ebp, msmo1, nsdhl, sqle, lss, and ogdh was significantly increased in the G20 group, while that of ebp, msmo1, and sqle was significantly decreased, compared with the G30 group. However, the mRNA expression of nsdhl, lss, and ogdh showed no

Validation of RNA Sequencing Data by QPCR
Considering the top 20 most important salinity stress responses of DEGs identified in the transcriptomic analysis, which had been reported in the previous studies, nine genes were selected for RNA sequencing data validation using qPCR ( Figure 6). The mRNA expression of ebp, msmo1, nsdhl, sqle, lss, and ogdh was significantly increased in the G20 group, while that of ebp, msmo1, and sqle was significantly decreased, compared with the G30 group. However, the mRNA expression of nsdhl, lss, and ogdh showed no difference between the G40 and G30 groups. The results were consistent with the RNA sequencing analysis results conducted using qPCR. Moreover, the mRNA expression of edar, wnt4, slc25a48 were decreased in the G20 group. However, the mRNA expression of edar and slc25a48 showed no difference between the G40 and G30 groups. Interestingly, the expression of these genes showed no differences at 15 days ( Figure S1).
with the G30 group. However, the mRNA expression of nsdhl, lss, and ogdh showed no difference between the G40 and G30 groups. The results were consistent with the RNA sequencing analysis results conducted using qPCR. Moreover, the mRNA expression of edar, wnt4, slc25a48 were decreased in the G20 group. However, the mRNA expression of edar and slc25a48 showed no difference between the G40 and G30 groups. Interestingly, the expression of these genes showed no differences at 15 days ( Figure S1). The relative expression level of mRNA transcripts was detected using qPCR via the 2 -∆∆Ct method. Data are expressed as means ± SE (n = 3). The asterisks * and ** indicate statistical differences at p < 0.05 and p < 0.01, respectively, as determined by one-way ANOVA with Tukey's post hoc test. βactin was used as the reference gene.

Transcriptomic Analysis of Differentially Expressed Genes
Salinity alterations can lead to various physiological reactions to maintain homeostasis, including osmotic regulation, ion transports, and respiratory metabolism [55]. For instance, low (12 ppt) and high salinity water (32 ppt) can significantly induce some specific pathways in gills of silvery pomfret (Pampus argenteus), including calcium transport, neuroactive ligand-receptor interaction, NOD-like receptor signaling, Toll-like receptor signaling, and cytokine-cytokine receptor interaction pathway, indicating that salinity stress affects the immune system and osmotic pressure-regulated pathways [30]. Moreover, low salinity stress (6 ppt) affects ion transport, immune response, energy metabolism, and protein synthesis in marbled flounder (Pseudopleuronectes yokohamae) [9]. A total of 417 DEGs (205 up-regulated and 212 down-regulated) in the 30 ppt vs. 20 ppt group and 18 DEGs (nine up-regulated and nine down-regulated) in the 30 ppt vs. 40 ppt group were identified. KEGG analysis showed that the DEGs were mainly enriched in the cytokinecytokine receptor interaction, apoptosis, steroid biosynthesis, and mTOR signaling path- The relative expression level of mRNA transcripts was detected using qPCR via the 2 −∆∆Ct method. Data are expressed as means ± SE (n = 3). The asterisks * and ** indicate statistical differences at p < 0.05 and p < 0.01, respectively, as determined by one-way ANOVA with Tukey's post hoc test. β-actin was used as the reference gene.

Transcriptomic Analysis of Differentially Expressed Genes
Salinity alterations can lead to various physiological reactions to maintain homeostasis, including osmotic regulation, ion transports, and respiratory metabolism [55]. For instance, low (12 ppt) and high salinity water (32 ppt) can significantly induce some specific pathways in gills of silvery pomfret (Pampus argenteus), including calcium transport, neuroactive ligand-receptor interaction, NOD-like receptor signaling, Toll-like receptor signaling, and cytokine-cytokine receptor interaction pathway, indicating that salinity stress affects the immune system and osmotic pressure-regulated pathways [30]. Moreover, low salinity stress (6 ppt) affects ion transport, immune response, energy metabolism, and protein synthesis in marbled flounder (Pseudopleuronectes yokohamae) [9]. A total of 417 DEGs (205 up-regulated and 212 down-regulated) in the 30 ppt vs. 20 ppt group and 18 DEGs (nine up-regulated and nine down-regulated) in the 30 ppt vs. 40 ppt group were identified. KEGG analysis showed that the DEGs were mainly enriched in the cytokine-cytokine receptor interaction, apoptosis, steroid biosynthesis, and mTOR signaling pathways, similar to cobia (Rachycentron canadum) [31]. As discussed below, the DEGs were involved in several potential complex molecular biological processes in the gill.

DEGs Involved in Steroid Biosynthesis and Lipid Metabolism
Steroid hormones, such as epinephrine and cortisol, can influence the metabolic capacity of the gill [56]. Epinephrine induces glycogenolysis after exposure to stress, thus increasing the plasma glucose level, which provides energy for the target tissue, including gill, to transfer ions like Na + [57]. In addition, cortisol can stimulate active Ca 2+ uptake under asymmetrical conditions, thus regulating the tight junction morphology between pavement cells of euryhaline fish [56,58]. Herein, KEGG pathway analysis showed that some DEGs (ebp, lss, sqle, nsdhl, msmo1, sc5d, dhcr7, and dhcr24) were involved in the steroid biosynthesis pathways under the long-term hyper-salinity (40 ppt) and hypo-salinity (20 ppt) stresses. Moreover, salinity stress (from 20 to 40 ppt) decreased the mRNA expression levels of sqle, msmo1, and ebp. Ebp (Emopamil binding protein), also known as EBP cholestenol delta-isomerase, is essential in the sterol biosynthesis pathway [59]. Sterols are essential cell membrane components and transporters in many biofilms [60]. SQLE (squalene epoxidase) is rate limiting and the first oxygenation enzyme in cholesterol synthesis [61]. Lanosterol, especially cholesterol, is the upstream precursor of sterol biosynthesis in fungal steroids and animals [62]. Herein, hypo-salinity up-regulated sqle, ebp, and lss in the gills, indicating the stimulation of cholesterol synthesis in the gills. Studies have shown that lss and dhcr24 are up-regulated in O. niloticus and Rachycentron canadum under hypo-salinity [31,63], consistent with this study.
Msmo1 is a key cholesterol biosynthetic enzyme. Nsdhl regulates adipogenesis via a synergized expression pattern with Msmo1 [64][65][66][67]. Herein, hypo-salinity up-regulated both nsdhl and msmo1, indicating that hypo-salinity can affect adipogenesis in greater amberjack. Moreover, hypo-salinity down-regulated edar and tnfsf12 mRNA in greater amberjack. Studies have shown that EDAr (ectodysplasin A receptor), belonging to the tumor necrosis factor receptor (TNFr) superfamily, can regulate cell activities, such as differentiation, proliferation, maturation, and lipid metabolism, by binding to the ectodysplasin A1 (EDA1) [68][69][70][71]. A previous study showed that tnfsf12 (TNF superfamily member 12) can inhibit lipid deposition in a dose-dependent manner without any cytotoxic effects. However, an agonistic antibody of the tnfsf12 receptor can alleviate the repression [72]. Herein, hypo-salinity decreased the mRNA expression of tnfsf12 in gill, indicating that lipid deposition is essential under low salinity stress. Previous reports have suggested that lipids are the energy source for euryhaline fish under osmotic stress [49,73]. Therefore, more lipids may be produced in greater amberjack under salinity changes by regulating the steroid biosynthesis related to lipid metabolism, adipogenesis, and lipid deposition, as revealed in other fish species [18]. The Ogdh enzyme is a key entry point for carbon into the Krebs cycle. Moreover, it affects all redox signals in mitochondria and cells [74,75]. Herein, hypo-salinity significantly increased the mRNA expression of ogdh mR, implying that adequate energy is needed in gills under salinity changes.

DEGs Involved in Ion Transport
Solute transport protein (solute carrier SLC) is the largest class of intracellular transport proteins, with over 300 members, mostly located in cell membranes. They mainly facilitate the transport of various substrates, including amino acids, nucleotides, glucose, and inorganic ions, across biological membranes [39]. Some of these proteins (slc5a6a, slc4a1a, slc4a4a) had different expression patterns between the hypo-salinity and control groups. The SLC4 gene family mediates HCO3 − extrusion and Cl − uptake across cellular plasma membranes, thus regulating the cell volume and intracellular pH and stabilizing resting membrane potential through the regulation of cytoplasmic Cl − [76]. Moreover, slc4a1a and slc4a1b are positively correlated with NKA (Na + /K + -ATPase) [77]. Herein, slc4a1a was down-regulated in the hypo-salinity group, while it was up-regulated in the hyper-salinity group, indicating the ion exchange in the gills.
Cystic fibrosis transmembrane conductance regulator (CFTR), an ABC transporter, acts as a channel across the cell membrane for transporting Cl − into and out of cells. Hypersalinity increases cftr expression in striped bass (Morone saxatilis) [78]. Moreover, cftr mRNA levels are significantly increased in gills of killifish (Fundulus heteroclitus) and Atlantic salmon (Salmo salar) under seawater conditions [79,80]. Similarly, this study showed that cftr was up-regulated in the hyper-salinity group, indicating the Cl − transport changes.
Aquaporins (AQPs) are a family of integral membrane proteins that facilitate water transport across biological membranes along an osmotic gradient [81]. Thirteen AQP isoforms (AQP0-AQP12) have been identified in humans and rodents [40]. However, AQP4 has the most potential for high water permeability [82]. Herein, the expression of aqp4 was significantly decreased in the hypo-salinity group compared with the control group, indicating the possible water permeability changes.

DEGs Involved in Immune Response
Previous studies have shown that the changes in fish immune status depend on the intensity and duration of the environmental stresses [83]. Herein, several DEGs (wnt4, slc25a48, slc6a8, tlr5, etc.) related to immune responses were found in the gill transcriptome of greater amberjack. Studies have reported that wnt genes enhance defense against pathogenic virus infection in the innate immune of Litopenaeus vannamei [84]. Moreover, Wnt4 protein can stimulate white blood cells and thymopoiesis in mice [85]. Herein, wnt4 was significantly down-regulated in the hypo-salinity group, indicating that low salinity can affect the immune system of the greater amberjack.
Toll-like receptors (TLRs) are key pathogen pattern recognition receptors that control the host immune responses against pathogens by recognizing molecular patterns specific to microorganisms [86]. Previous studies have reported that tlr5 plays a crucial role in the immune responses of turbot (Scophthalmus maximus L.) to the infections of various pathogens [87]. For instance, TLR5 stimulates the expression of proinflammatory, antibacterial, and stress-related genes by binding to bacterial flagellin, thus enhancing host defense against bacterial pathogens [88]. Moreover, flagellin increases Tlr5 activity and the release of its downstream factor, IL−8 in mice [88,89]. Herein, the expression level of tlr5 was decreased in the hypo-salinity group, indicating that low salinity inhibits tlr5 activity on key adaptive functions, thus lowering efficient immune responses.

Conclusions
This work used transcriptome analysis to investigate the molecular changes in gills of greater amberjack (Seriola dumerili) under three different salinity concentrations (20, 30, and 40 ppt). The results provided large transcriptome data, abundant DEGs, and signaling pathways related to salinity adaptation. The signaling pathways analysis indicated that a complex molecular regulatory network is involved in metabolism, including steroid synthesis, lipid metabolism, tryptophan metabolism, ion transporters, and immune response for adaptation to salinity stress. Therefore, this study can provide insights into the molecular mechanisms of greater amberjack adaptation to salinity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani12030327/s1, Figure S1: Relative expression levels of different genes in RT-PCR (samples taken on the 15th day); Table S1: Primer sequences used in this study.

Conflicts of Interest:
The authors declare no conflict of interest.