Abstract

   Acanthogobius ommaturus is a euryhaline fish widely distributed in
   coastal, bay and estuarine areas, showing a strong tolerance to
   salinity. In order to understand the mechanism of adaptation to
   salinity stress, RNA-seq was used to compare the transcriptome
   responses of Acanthogobius ommaturus to the changes of salinity. Four
   salinity gradients, 0 psu, 15 psu (control), 30 psu and 45 psu were set
   to conduct the experiment. In total, 131,225 unigenes were obtained
   from the gill tissue of A. ommaturus using the Illumina HiSeq 2000
   platform (San Diego, USA). Compared with the gene expression profile of
   the control group, 572 differentially expressed genes (DEGs) were
   screened, with 150 at 0 psu, 170 at 30 psu, and 252 at 45 psu.
   Additionally, among these DEGs, Gene Ontology (GO) analysis indicated
   that binding, metabolic processes and cellular processes were
   significantly enriched. Kyoto Encyclopedia of Genes and Genomes (KEGG)
   pathways analysis detected 3, 5 and 8 pathways related to signal
   transduction, metabolism, digestive and endocrine systems at 0 psu, 30
   psu and 45 psu, respectively. Based on GO enrichment analysis and
   manual literature searches, the results of the present study indicated
   that A. ommaturus mainly responded to energy metabolism, ion transport
   and signal transduction to resist the damage caused by salinity stress.
   Eight DEGs were randomly selected for further validation by
   quantitative real-time PCR (qRT-PCR) and the results were consistent
   with the RNA-seq data.

   Keywords: Acanthogobius ommaturus, RNA-seq, salinity stress,
   differentially expressed genes (DEGs)

1. Introduction

   The salinity is one of the most significant factors for regulating
   distribution, abundance and diversity of aquatic animals [[30]1], and
   it is also one of the main environmental factors that exerts selective
   pressure on aquatic organisms. The sea surface salinity would rise with
   the increase of temperature under global warming [[31]2] Changes in
   environmental salinity will directly affect osmotic pressure of aquatic
   organisms [[32]3], which has been predicted to have a profound impact
   on fish at the molecular and cellular level. Fish need to make adaptive
   changes to maintain physiological function, and conduct compensatory
   adjustment through changing distribution range [[33]4].

   Acanthogobius ommaturus is a large, demersal fish of the Gobiidae,
   which is economically important across its distribution area [[34]5].
   It is a euryhaline fish that could adapt rapidly and maintain
   homeostasis in a wide range of salinities (from freshwater to salinity
   30) [[35]6]. Studies on A. ommaturus mainly focused on the comparative
   analysis of individual fecundity and morphology [[36]7], fishery
   biology [[37]8], name of species [[38]9], and genetic diversity and
   structure by molecular markers [[39]10]. However, the potential
   mechanism of salinity adaptation in A. ommaturus, especially at the
   molecular level, such as comprehensive molecular pathway response, has
   not been well studied. It is worth noting that the Acanthogobius
   flavimanus, which is closely related to the A. ommaturus and
   distributed in East Asia, has invaded California, USA [[40]11,[41]12]
   and southern Australia [[42]13]. As a euryhaline fish with strong
   viability and rapid growth, A. ommaturus may be also capable to invade
   other waters leading to the imbalance or collapse of the local marine
   ecosystem. Therefore, it is very necessary to conduct the potential
   mechanism of salinity adaptation of A. ommaturus.

   With the development of high-throughput sequencing technology, RNA-seq
   provides an opportunity to study transcriptome variation in organisms
   under different environmental conditions. Transcriptome sequencing can
   not only detect almost all effective genes expressed in specific cells
   or organs, but also provide a comprehensive understanding of the
   regulatory mechanisms involved in specific biological processes based
   on the structure and function of differential genes [[43]14]. In
   addition, transcriptome sequencing allows simultaneous analysis of all
   physiological process, including metabolism [[44]15], protein
   homeostasis [[45]16] and other regulatory cellular processes [[46]17].
   Thus, transcriptome technology addresses key gene expression patterns
   of non-model organisms in specific environments and enables detection
   of unknown genes and discovery of new transcripts. In recent years,
   transcriptome technology has been widely used in transcriptome assembly
   and annotation of many non-model fish, such as Megalobrama amblycephala
   [[47]18], Pelteobagrus fulvidraco [[48]19] and Atlantic salmon
   [[49]20].

   In this study, to further know about the osmoregulatory mechanisms of
   A. ommaturus in response to salinity stress, RNA-Seq was used to
   sequence transcripts of gill tissues of A. ommaturus under salinity
   stress. This is the first transcriptome study of A. ommaturus, which
   will provide important information for enriching its genetic resources.
   Moreover, the main purpose of this study is to explore and identify the
   genes and pathways of gill organs that play a potential role in salt
   tolerance by using RNA-Seq.

2. Materials and Methods

2.1. Ethics Statement

   All animal experiments were conducted in accordance with the guidelines
   and approval of the respective Animal Research and Ethics Committees of
   Ocean University of China. The field studies did not involve endangered
   or protected species. In addition, frost anesthesia was performed to
   minimize the suffering of all A. ommaturus specimens.

2.2. Experimental Design and Sampling

   Acanthogobius ommaturus specimens (initial weight 41.88 ± 5.52 g,
   initial length 19.23 ± 1.91 cm) were collected from the coastal water
   of Qingdao, China. All specimens were dispatched into four aquariums
   (80 × 60 × 40 cm) with recirculating aerated natural seawater (25 psu)
   at a density of 20 fishes per aquarium for one week of acclimation.
   During the acclimation and experimental periods, A. ommaturus was fed
   to satiation with small shrimp twice daily (08:30 and 16:00), and the
   remaining feed and feces were siphoned out and 1/3 of the water in each
   aquarium was replaced with pre-aerated water each day. After 7 days of
   acclimation, healthy fish with strong individual vitality were used for
   salinity stress experiments. In the manner of acute stress, each box
   was a salinity gradient, and the four salinity gradients were 0 psu, 15
   psu, 30 psu and 45 psu, respectively. In this experiment, since the
   optimum salinity range of A. ommaturus is 5–15 psu [[50]21], salinity
   15 was selected as the control salinity. Three groups of experiments
   with salinity of 0 vs. 15, 15 vs. 30 and 15 vs. 45 were analyzed. After
   48 h salinity stress, three experimental individuals with similar size,
   vigor and health were immediately euthanized and sampled in each
   experimental group. Gills tissues were collected and stored at −80 °C
   for RNA extraction.

2.3. RNA Isolation and Illumina Sequencing

   Total RNA was extracted from the gill tissues using the TRIzol reagent
   (Invitrogen, Carlsbad, USA) according to the instructions, and genomic
   DNA was removed using DNase I (TaKaRa, Otsu, Japan). RNA purity was
   checked using the NanoPhotometer^® spectrophotometer (IMPLEN,
   Calabasas, CA, USA). RNA concentration was measured using Qubit^® RNA
   Assay Kit in Qubit^® 2.0 Flurometer (Life Technologies, Carlsbad, CA,
   USA). RNA quality was determined by an Agilent 2100 Bioanalyzer. Only
   high-quality RNA samples (OD260/280 = 1.9–2.1, OD260/230 ≥ 2.1, RIN ≥
   9.5, 28/18S ≥ 1.0, > 30 μg) were used to construct these quenching
   libraries. The RNAs from tissues of salinity gradient were pooled in
   equal amounts and mRNA was extracted from the total RNA using magnetic
   beads with Oligo (dT) probes. Fragmentation buffer was applied to lyse
   the mRNA into fragments with a suitable size and the fragmented mRNA
   was used to construct a cDNA library using TruSeq Stranded mRNA
   LTSample Prep Kit (Illumina, San Diego, CA, USA). Then, the sequencing
   was carried out by Novogene Company (Beijing, China) on Illumina Hiseq
   2000 platform.

2.4. De Novo Assembly and Functional Annotation

   All raw reads in FASTQ format quality control relied on NGS QC ToolKit
   software (v2.3.3, National Institute of Plant Genome Research, New
   Delhi, India) [[51]22] and clean reads were obtained by removing raw
   reads containing adapter, ploy-N (N ratio > 10%) and low quality reads
   (quality scores ≤ 5). The Q20, Q30, GC (Guanine and Cytosine)-content
   and sequence duplication level of the clean data were calculated. All
   the downstream analyses were based on clean data with high quality.
   Transcriptome de novo assembly was carried out with Trinity software
   (version 2.0.6 Broad Institute, Cambridge, MA USA) with a min_kmer_cov
   set to 2 by default and all other parameters set to their defaults
   [[52]23]. After assembly, the transcripts were clustered and the
   redundant ones were removed and the remaining sequences were defined as
   unigenes. All of the assembled unigenes were searched against the Nr
   (NCBI non-redundant protein sequences) database, the eukaryotic
   orthologous group (KOG) database, COG (Clusters of Orthologous Groups
   of Proteins) database, Pfam (Protein family) database, SwissProt
   databases, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database
   and Gene Ontology (GO) database using BLASTx to identify the proteins
   that had the highest sequence similarity with the given transcripts to
   retrieve their function annotations, and a typical E-value cut-off was
   set at < 1.0 × 10^−5 [[53]24].

2.5. Differential Expression Analysis

   To investigate the different gill transcriptome responses of A.
   ommaturus to salinity fluctuation, we analyzed the number and
   biological functions of differentially expressed genes (DEGs) in three
   treatment groups. Firstly, we used BWA-mem [[54]25] to map all unigenes
   to a multicast file. Then, the expression levels of all unigenes were
   normalized using RSEM [[55]26] and Bowtie2 [[56]27] to determine the
   number of fragments per kilobase of exon model per million mapped
   fragments (FPKM). DEGs of three treatment salinity were identified
   using the edgeR package [[57]28], and |log[2]FC| ≥ 1 (which corresponds
   to FC = 2) and FDR ≤ 0.05 were used as the filtering thresholds. In
   addition, in order to further understand the function of DEGs under
   salinity stress, DEGs were used to determine the Gene Ontology (GO)
   term and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.
   Blast2GO ([58]https://www.blast2go.com/) software was used to obtain
   the GO annotations of unigenes for describing biological processes,
   molecular functions and cellular components [[59]29]. The KEGG
   Automatic Annotation Server (KAAS)
   ([60]http://www.genome.jp/kaas-bin/kaas_main) system was used for
   pathway reconstruction, respectively [[61]30].

2.6. Quantitative Real-Time PCR

   In order to verify the reliability of the transcriptome sequencing
   results, eight differentially expressed genes (DEGs) with high
   expression levels were randomly selected, and primers were designed
   with Primer 6.0 ([62]Table S1). Firstly, reverse transcription was
   carried out from each total RNA sample using the PrimeScript RT reagent
   Kit with gDNA Eraser (TaKaRa, JPN) with RT Primer Mix of Random 6 mers
   and Oligo dT Primer. Furthermore, standard curves were constructed to
   identify the ideal dilution times of cDNA samples and were used as
   calibrators. A total of 12 cDNA samples were diluted in 10-fold with
   nuclease-free water and were used as templates for PCR. Furthermore,
   the qRT-PCR analysis was designed following instructions for SYBR^®
   Premix Ex TaqTM (Tli RNaseH Plus) RR420A with a StepOnePlus (TaKaRa,
   JPN). A reaction system of 20 μL was amplified, including 2.0 μL of
   diluted cDNA template, 10μL of SYBR Premix Ex Taq (2×), 0.4 μL of each
   of the forward and reverse primers and ROX Reference Dye (50×) and 6.8
   μL of nuclease-free water. The amplification processes consisted of a
   holding stage of 30 s at 95 °C, followed by 40 cycles of 5 s at 95°C
   and 35 s at 52 °C. All reactions were performed in triplicates. Melting
   curve analysis was performed to determine the target specificity, and
   melting curve temperature was 61.9–86.7 °C. The relative expression
   levels of all target unigenes were calculated by the 2^−ΔΔCT analysis
   method (ΔCT = CT[target unigene] − CT[reference gene], ΔΔCT =
   ΔCT[treatment] − ΔCT[control]), and β-actin was used as the reference
   gene for qRT-PCR normalization.

2.7. Statistical Analysis

   Statistical analyses were performed using SPSS 18.0 software. Data
   expressed as the mean ± standard deviation (SD) from three independent
   replicates. Significant differences between samples were analyzed by
   Duncan’s tests at a significance level of 0.05.

3. Results

3.1. Illumina Sequencing and Assembly

   In brief, RNA-Seq was carried out on libraries constructed from gill of
   A. ommaturus. After quality trimming and adapter clipping, results from
   RNA-Seq produced 12 libraries that consisted of 638,902,551 clean reads
   with an N50 value equal to 1111. The error rate was less than 0.03,
   Q30% (the rate of bases for which quality is greater than 20) was over
   93%, and GC percentage was around 47–49%, indicating high-quality
   sequencing and the feasibility of the subsequent analyses. After de
   novo assembly using trinity software (version 2.0.6 Broad Institute,
   Cambridge, MA, USA), 131,225 unigenes and 88,228,207 assembled bases
   were assembled. The overall sequencing results are shown in [63]Table
   S2. All raw sequence data were deposited in NCBI Sequence Read Archive
   (SRA) under accession number SRP258637.

3.2. Unigenes Functional Annotation and Classification

   Unigenes were subjected to annotation analysis by comparing with Nr,
   KOG, COG, SwissProt, Pfam, KEGG, and GO databases. Results show that a
   total of 46,443 unigenes (35.38%) were annotated in at least one
   database, with 39,920 annotated unigenes (30.42%) had a significant
   BLAST hit against Nr database. The detailed annotation results are
   listed in [64]Table 1. For top-hit species matched against Nr database,
   51.4% of the matched unigenes indicated similarity with Boleophthalmus
   pectinirostris (51.4%), followed by Larimichthys crocea (3.8%) and
   others (31.9%) ([65]Figure S1).

Table 1.

   Summary of annotation of A. ommaturus unigenes.
   Database Number of Annotated Unigenes Percentage of Annotated Unigenes
   Nr 39,920 30.42%
   KOG 21,947 16.72%
   COG 11,807 9.00%
   Pfam 27,697 21.11%
   SwissProt 20,465 15.60%
   KEGG 22,777 17.36%
   GO 18,924 14.42%
   At least one database 46,443 35.38%
   Total unigenes 131,225 100.00%
   [66]Open in a new tab

3.3. Analysis of Differentially Expressed Genes (DEGs)

   To elucidate the gene expression pattern under salinity fluctuation, we
   firstly compared the numbers of DEGs of three treatment groups. A total
   of 150 (59 up- and 91 down-regulated), 170 (39 up- and 131
   down-regulated), and 252 (159 up- and 93 down-regulated) genes were
   differently expressed with |log[2]FC|≥1 (which corresponds to FC = 2)
   and FDR ≤ 0.05 in 0 vs. 15 psu, 15 vs. 30 psu and 15 vs. 45 psu,
   respectively. As shown in [67]Figure 1, three DEGs were present in all
   three comparisons.

Figure 1.

   [68]Figure 1
   [69]Open in a new tab

   Global overview of differentially expressed unigenes based on Venn
   diagram.

3.4. GO and KEGG Pathway Enrichment Analysis of DEGs

   To identify the functional changes potentially associated with the
   salinity changes of A. ommaturus, GO and KEGG pathway enrichment
   analyses were conducted. Gene ontology analysis of the three
   comparisons was performed to ascribe GO categories to 150, 170 and 252
   DEGs, respectively. Compared to the control group (15 psu), 101 of 150,
   109 of 170 and 167 of 252 significant DEGs were classified into 37 GO
   terms, 34 GO terms and 39 GO terms at 0 psu, 30 psu and 45 psu,
   respectively. The distribution of GO terms indicated that cellular
   processes (GO: 0009987) and metabolic processes (GO: 0008152) in
   biological processes, cells (GO: 0005623) and cell parts (GO: 0033643)
   in cellular components, and binding (GO: 0005488) and catalytic
   activity (GO: 0003824) in molecular functions were most significantly
   enriched under salinity stress. The detailed information regarding
   significant GO terms is shown in [70]Figure 2.

Figure 2.

   [71]Figure 2

   [72]Figure 2
   [73]Open in a new tab

   Results of the enrichment analysis of gene ontology (GO) terms of
   differentially expressed genes (DEGs) of A. ommaturus exposed to
   salinity changes.

   KEGG pathways provide valuable information for studying the specific
   biological, metabolic processes and molecular mechanisms under salinity
   stress in A. ommaturus. With the increase of salinity, DEGs were
   enriched in 170, 168 and 234 different pathways, among which 7, 9 and
   11 pathways (q value < 0.05) were significantly enriched in 0 vs. 15,
   15 vs. 30 and 15 vs. 45 groups, respectively. Compared to the control
   group, most DEGs were enriched in signal transduction pathways,
   metabolic pathways, digestive and endocrine systems at 0 psu, 30 psu
   and 45 psu, respectively. For the KEGG pathway analysis, the dominant
   pathways were as follows: cytokine–cytokine receptor interaction
   (ko04060), metabolism of xenobiotics by cytochrome P450 (ko00980) and
   proximal tubule bicarbonate reclamation” (ko04964) in 0 vs. 15, 15 vs.
   30 and 15 vs. 45, respectively. The significantly enriched pathways in
   the top 20 were shown in [74]Figure 3.

Figure 3.

   [75]Figure 3

   [76]Figure 3
   [77]Open in a new tab

   Scatterplot of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways
   enriched in the DEGs. Rich factor is the ratio of the DEGs number to
   the total gene number in a given pathway. The size and color of the
   dots represent the gene number and range of the q value, respectively.

3.5. Identification of DEGs Related to Salinity Changes

   As a non-model species with only limited gene function annotation
   resources, we based Nr annotation and manual literature searches on the
   combination of enrichment analysis to further discuss candidate DEGs
   potentially associated with salinity adaptation and osmoregulation.
   These DEGs were categorized into three functional categories including
   energy metabolism, ion transporters and signal transduction ([78]Table
   2) by GO enrichment analysis. Imputed putative functions of these genes
   are covered in the discussion.

Table 2.

   Representative osmotic pressure related genes involved in salinity
   stress.
   Functional Group Group Gene Name Gene Function Gene ID log[2] (Fold
   Change)
   Energy metabolism 0 vs. 15 glutathione S-transferase(rho)(GST)
   transferase activity 96847_c4_g1 −1.117
   mevalonate kinase isoprenoid biosynthetic process 94143_c2_g1 1.148
   insulin-like growth factor-binding protein 1(IGFBP1) regulation of cell
   growth 78537_c0_g1 −1.489
   inositol-tetrakisphosphate 1-kinase-like inositol trisphosphate
   metabolic process 94029_c1_g1 1.982
   15 vs. 30 Calreticulin (CRT) protein folding 88942_c1_g2 −1.136
   transmembrane protease serine 9 proteolysis 95486_c1_g3 −1.234
   glutamine synthetase (GS) glutamine biosynthetic process 85559_c6_g7
   1.754
   phosphoenolpyruvate carboxykinase(PEPCK) gluconeogenesis 89938_c3_g1
   1.763
   glutathione S-transferase(rho)(GST) transferase activity 96847_c4_g1
   -1.907
   glutamate decarboxylase 1(GAD1) carboxylic acid metabolic process
   94253_c5_g1 -1.289
   15 vs. 45 glutamine synthetase (GS) glutamine biosynthetic process
   65875_c1_g1 1.141
   serine/threonine-protein kinase Sgk1 signal transduction mechanisms
   87076_c4_g1 1.347
   Ion transporters 0 vs. 15 sodium potassium ATPase beta (NKA) potassium
   and sodium ion transport 92532_c3_g1 1.245
   transient receptor potential cation channel subfamily V member 4(TRPV4)
   ion transmembrane transport 90898_c2_g1 1.341
   Aquaporin-3(AQP3) water transport 94106_c4_g1 2.08
   FXYD domain-containing ion transport regulator 11(FXYD11) ion transport
   97290_c7_g1 −1.606
   sodium-coupled neutral amino acid transporter 2(SNAT2) amino acid
   transport and metabolism 91967_c1_g2 −1.087
   cystic fibrosis transmembrane conductance regulator (CFTR) chloride
   transmembrane transport 89496_c1_g1 −1.458
   Na^+:K^+:2Cl^− cotransporter 1a(NKCC1a) Na^+, K^+, and Cl^- transport
   92335_c1_g1 −1.004
   15 vs. 30 ferritin middle subunit iron ion transport 90148_c4_g2 −1.265
   carbonic anhydrase 4(CA4) metal ion binding 96156_c2_g1 1.480
   adenosylhomocysteinase 3 isoform X7 sodium ion transport 96455_c7_g1
   1.007
   hemoglobin subunit α-1 (HBA1) oxygen transport 84032_c0_g1 −1.764
   aquaporin-1 (AQP1) water transport 82312_c0_g1 −2.372
   15 vs. 45 Na^+/K^+ ATPase α1b-ii (NKA) potassium and sodium ion
   transport 93985_c4_g5 1.745
   sodium-coupled neutral amino acid transporter (SNAT) amino acid
   transport and metabolism 91115_c3_g4 1.893
   aquaporin 3(AQP3) water transport 94106_c4_g1 1.532
   FXYD domain-containing ion transport regulator 11(FXYD11) ion transport
   97290_c7_g1 1.791
   chloride channel protein 2(CLCP2) chloride transmembrane transport
   79986_c0_g1 −1.134
   Na^+:K^+:2Cl^− cotransporter 1(NKCC1a) Na^+, K^+, and Cl^- transport
   92335_c1_g1 1.495
   cystic fibrosis transmembrane conductance regulator (CFTR) chloride
   transmembrane transport 89496_c0_g1 1.202
   FXYD domain-containing ion transport regulator 9(FXYD9) ion
   transmembrane transport 95234_c8_g2 1.018
   sodium/potassium-transporting ATPase potassium ion transport
   85132_c6_g1 1.762
   ATP-sensitive inward rectifier potassium channel 15(KCNJ15) potassium
   ion import 68888_c0_g1 1.829
   Signal transduction 0 vs. 15 phosphatidylinositol 3-kinase regulatory
   signal transduction 79932_c2_g2 1.092
   phospholipid phosphatase 3-like isoform X1 signal transduction
   93118_c0_g1 1.062
   cAMP-responsive element modulator (CREM) regulation of transcription
   85301_c4_g2 −2.242
   cAMP-specific 3’,5’-cyclic phosphodiesterase 4B(PDE4B) signal
   transduction 88103_c7_g1 −1.041
   prolactin receptor (PRLR) prolactin signaling pathway 92757_c8_g1 2.672
   15 vs. 30 high affinity cAMP-specific and IBMX-insensitive 3’,5’-cyclic
   phosphodiesterase 8A signal transduction 83321_c0_g1 −1.073
   putative transmembrane protein 116(PTP) G-protein coupled receptor
   signaling pathway 85145_c0_g1 −2.122
   15 vs. 45 GTPase HRas small GTPase mediated signal transduction
   91669_c2_g4 1.139
   cAMP-responsive element modulator (CREM) regulation of transcription
   85301_c4_g2 1.089
   cyclooxygenase-2 (COX2) inflammatory response 78208_c1_g1 1.132
   cAMP-specific 3’,5’-cyclic phosphodiesterase 4B(PDE4B) signal
   transduction 90244_c2_g1 1.445
   Inositol monophosphatase(IMP) signal transduction 94682_c12_g1 1.652
   MAP kinase-interacting serine/threonine-protein kinase 2-like
   intracellular signal transduction 91691_c2_g1 1.333
   [79]Open in a new tab

3.6. Validation of Transcriptomic Data by qRT-PCR

   Eight differentially expressed genes (DEGs) were randomly selected to
   validate the reliability of the transcriptome in this study.
   Melting-curve analysis revealed a single product for all of the tested
   genes. For these candidate genes, the variation trend in expression was
   concordant between qRT-PCR data and transcriptomic data, excluding two
   genes (86284_c2_g1 in 0 vs. 15 group and 91346_c0_g1 in 15 vs. 30
   group) did not match perfectly. Therefore, the results indicated the
   transcriptomic data were credible ([80]Figure 4).

Figure 4.

   [81]Figure 4
   [82]Open in a new tab

   Relative change of the transcriptomes data and qRT-PCR data of 8 DEGs.
   Log[2]fold changes are expressed as the ratio of gene expression after
   normalization to β-actin. The data are means ± SD from three
   independent replicates, and “*” indicates statistical significance
   between experimental group and control group (p < 0.05).

4. Discussion

   Salinity is one of the decisive environmental factors affecting fish
   survival and growth [[83]31], food intake [[84]32], disease [[85]33]
   and distribution [[86]34]. However, limited evidence is available to
   explain how salinity stress affects the cellular physiological effects
   of A. ommaturus. RNA-Seq analysis is considered to be a robust method
   for evaluating transcriptional responses to different experimental
   conditions, especially in non-model organisms without reported genomes
   [[87]35]. In this study, 156, 165, 167 and 151 million clean reads were
   obtained for four salinity gradients of 0 psu, 15 psu, 30 psu and 45
   psu, respectively. For most eukaryotic transcriptomes, sequencing up to
   100 million reads can quantify precisely genes and transcripts that
   have low expression levels [[88]36]. According to the analysis of RNA
   sequencing data, Q30 is more than 93%, and the N50 length of unigenes
   is 1111 bp, which proves the reliability of raw data. From the
   distribution of species annotated by nr database, 20,505 of 39,920
   unigenes (51.4%) have been annotated to the species of Boleophthalmus
   pectinirostris. This is probably due to the closer phylogenetic
   relationship between two species since both are in the order of
   Gobioidei, which is also supported the accuracy of raw data. In
   addition, to further confirm gene expression profiles, qRT-PCR analysis
   of eight selected DEGs was conducted. Although there were some
   differences in gene expression between qRT-PCR and RNA-Seq, the overall
   gene expression trends were consistent. In conclusion, the sequencing
   results of this study are accurate and reliable, which is suitable for
   further analysis.

   To further explore the research results, we performed GO and KEGG
   enrichment analysis of differential expressed genes (DEGs). The results
   of GO enrichment analysis indicated that DEGs were significantly
   enriched in binding (GO: 0005488), metabolic processes (GO: 0008152),
   and catalytic activities (GO: 0003824). Previous studies have confirmed
   these three functions play an important role in osmoregulatory
   mechanisms. Firstly, the 14-3-3 protein inhibit Ca^2+ activated Cl^−
   channels by binding (GO: 0005488) to calmodulin to maintain osmotic
   pressure homeostasis in the euryhaline teleost Fundulus heteroclitus
   [[89]37]. Acipenser brevirostrum appeared to have successfully adapted
   to higher salinity, probably due to the osmotic pressure balance
   maintained by the metabolism (GO: 0008152) of non-esterified fatty
   acids [[90]38]. In addition, previous studies had shown that the
   catalytic function (GO: 0003824) of enzymes regulated ion changes and
   osmotic pressure when crustaceans were exposed to salinity fluctuations
   [[91]39]. KEGG pathway analysis can enhance our understanding of how
   salinity fluctuations affect the osmoregulation of A. ommaturus. For
   the 0 vs. 15 group, cytokine receptor interaction and JAK STAT
   signaling pathway related to signal transduction was upstream pathways
   regulating immune and osmoregulation [[92]40]. Steroid hormone
   biosynthesis and amino acid metabolism associated with metabolic
   pathways involved in osmoregulation were significantly enriched in the
   15 vs. 30 group. Steroid hormone biosynthesis in metabolic pathways
   played pivotal roles in response to salinity stress and was
   inextricably linked to other pathways in aquatic animals [[93]41]. The
   variation of amino acids will directly lead to the change of protein
   content [[94]38], which facilitates more rapid acclimation to
   environmental changes than relying solely on enzymatic mechanisms
   [[95]42]. Regarding the annotated path map ([96]Figure S2) in the 15
   vs. 45 group, pancreatic secretion, mineral absorption, proximal tubule
   bicarbonate reclamation and bile secretion have a common feature of
   increasing the ATPase or ATP content to provide energy for the exchange
   of sodium and potassium ions. In these pathway maps, we found high
   expression of genes involved in osmotic pressure regulation, such as
   cystic fibrosis transmembrane conductance regulator (CFTR),
   Na^+-K^+-2Cl^- cotransporter 1 (NKCC1), and phosphoenolpyruvate
   carboxykinase (PEPCK). In addition, the energy metabolism, ion
   transporters and signal transduction related DEGs and their potential
   functions of salinity response were discussed below.

4.1. DEGs Related to Energy Metabolism

   In order to maintain osmotic balance, the osmoregulatory processes are
   regulated by several kinds of enzymes, and the synthesis and operation
   of these proteins requires large amounts of energy [[97]4]. Stress
   induced by salinity changes has been associated with enhanced reactive
   oxygen species (ROS) generation, causing oxidative damage [[98]43].
   Glutathione S-transferase (GST) is a key enzyme of cellular
   detoxification systems that defends cells against ROS [[99]44].
   However, compared with the control salinity, GST was down-regulated at
   0 psu and 30 psu in this study. These differences in GST response may
   highlight differences in detoxification capacity between species and
   different tissues [[100]45]. We hypothesized that the down-regulated
   expression of GST at 0 psu and 30 psu was related to the salinity
   tolerance of A. ommaturus. Glutamine synthetase (GS) activity may
   facilitate the production of glutamine, proline, and other organic
   solutes characteristic of osmotic and pH adjustments under stress
   conditions [[101]46]. GS expression was significantly up-regulated
   under high salinity stress (30 psu and 45 psu). We speculated that the
   high expression of GS promoted the production of proline and other
   amino acids, which was widely used in hyperosmotic fermentation systems
   as an osmotic pressure protectant [[102]47]. It is worth noting that
   phosphoenolpyruvate carboxykinase (PEPCK) was up-regulated at 30 psu.
   PEPCK catalyzes the conversion of oxaloacetate into
   phosphoenolpyruvate, which is the first rate-limiting step in
   gluconeogenesis [[103]48]. The increased expression of PEPCK in the
   gills of A. ommaturus indicated that the endogenous glucose increased
   by enhancing the gluconeogenesis to provide the extra energy to
   maintain the ion balance under the change of salinity [[104]48]. In
   brief, we hypothesized that genes associated with energy metabolism
   provided energy for cell survival and ultimately increased the ability
   of A. ommaturus to survive in salinity changes.

4.2. DEGs Related to Ion Transporters

   Osmotic balance is maintained by the coordinated transport between
   water and ions in the intestines, gills and kidneys in teleosts
   [[105]49]. The gill is an osmoregulartory organ; several classical ion
   transporters and channels were found in this study, such as Aquaporin
   (AQP1 and AQP3), transient receptor potential cation channel subfamily
   V member 4 (TRPV4), Na^+-K^+-2Cl^- cotransporter 1(NKCC1),
   ATP-sensitive inward rectifier potassium channel 15(KCNJ15) and sodium
   potassium ATPase beta subunit (NKA) [[106]50,[107]51,[108]52,[109]53].
   Tilapia is one of the typical experimental models for euryhaline fish.
   It would be more informative if we discuss the similarity and
   differences in the patterns of gene expression of gill ion transporters
   and channels between tilapia and A. ommaturus. To acclimate in a high
   salinity environment, NKA transports 3Na^+ outward in exchange for
   2K^+, creating low intracellular Na^+ and a highly negative charge
   within the cell; the Na^+ gradient is used to transport Na^+, K^+ and
   2Cl^− into the cell through a basolateral NKCC1; Cl^− then leaves the
   cells down an electrical gradient through an apical CFTR [[110]54]. In
   this osmoregulatory mechanism, the up-regulation of NKA, NKCC1a and
   CFTR in the gills of A. ommaturus does make sense; the similarity in
   the patterns of gene expression in Oreochromis mossambicus [[111]55]
   adapted to seawater stress. The currently accepted model for active
   NaCl secretion by ionocytes in seawater-acclimated teleosts consists
   primarily of the cooperative action of three major ion-transport
   proteins: NKA, NKCC1, and CFTR [[112]54]. Under low salinity
   environment, the basolateral NKA activity is probably involved in
   driving uptake of NaCl, possibly in conjunction with an apical V-type
   H^+-ATPase, via apical Na^+ channels (90406_c5_g1, log[2]FC = 0.231)
   and Cl^–/HCO[3]^– exchangers [[113]56]. The up-regulation of NKA plays
   an important role in freshwater adaptation of A. ommaturus and O.
   mossambicus [[114]55]. In the same way, the up-regulation of NKA has
   been found in Salvelinus namaycush and Salvelinus fontinalis [[115]57],
   and in Stenogobius hawaiiensis [[116]50] in both high and low salinity
   environments. NKA provides the major driving force for both ion
   secretion and absorption [[117]58], this may be closely related to the
   up-regulated expression of NKA in both high and low salinity
   environments.

   Teleost fishes are able to acclimatize to seawater by secreting excess
   NaCl by means of specialized “ionocytes” in the gill epithelium
   [[118]59]. The main function of NKCC is ion secretion
   [[119]59,[120]60], which explains that the up-regulation of NKCC under
   high salinity environment and down-regulated under low salinity
   environment in A. ommaturus. Unlike A. ommaturus, the up-regulation of
   NKCC in O. mossambicus occurred during freshwater adaptation, which was
   surprising [[121]55]. Moreover, ion absorption function Na^+/Cl^–
   cotransporter (NCC) and Na^+/H^+ exchanger 3 (NHE3), were up-regulated
   in freshwater adaptation of O. mossambicus [[122]58,[123]60], but this
   significant phenomenon was not found in A. ommaturus. Although the
   expression of NHE3 in A. ommaturus was up-regulated in the low salinity
   environment, the log[2]FC value was only 0.167. The similarity and
   differences in the patterns of gene expression in salinity adaptation
   between A. ommaturus and O. mossambicus suggests the likelihood that
   there is more than one type of ionocyte in most teleost fishes
   adaptable to freshwater, while ion-secretory ionocytes in seawater seem
   to be of only one type in all teleosts [[124]55]. Prolactin (PRL) is a
   freshwater adapting hormone in teleosts that acts by reducing water
   permeability and promoting ion uptake by osmoregulatory organs, such as
   gills [[125]61]. PRL release in a hypoosmotic environment is induced by
   an increase in PRL cell volume brought about by the influx of water via
   a membrane-bound AQP3 [[126]62] and Ca^2+ influx through the
   stretch-activated channel TRPV4 [[127]63]. PRL may play an important
   role in the control of water and electrolyte balance through PRLR
   expressed in the osmoregulatory organs in the marine teleost [[128]64].
   Despite, the significantly up-regulated expression of PRLR instead of
   PRL, we suspected that the up-regulated expression of AQP3 and TRPV4 at
   0 psu were the direct reason that promoted the secretion of PRL and led
   to the up-regulated expression of PRLR to adapt to freshwater stress.
   We found that the expression patterns of AQP3 and TRPV4 genes in A.
   ommaturus and O. mossambicus were similar, both of which were
   up-regulated under freshwater stress [[129]62,[130]63]. In addition,
   AQP1 in gills of Acanthopagrus schlegeli [[131]65] and Dicentrarchus
   labrax [[132]66] decreased with the increase in salinity, which was
   consistent with the results of this experiment. In the same way, FXYD
   also plays an important role in osmoregulation [[133]67]. In general,
   the patterns of gene expression of osmoregulatory mechanisms between A.
   ommaturus and O. mossambicus were largely similar. The synergistic
   effect of ion transportations and channels elicits complex regulatory
   strategies and ultimately increases tolerance to salinity stresses.

4.3. DEGs Related to Signal Transduction

   The adaptability and stimulatory response of fish to salt stress
   depends on the effective mechanism of osmosensing and osmotic stress
   signaling [[134]68]. Furthermore, osmosensory signals are not only
   rapidly transduced within cells, but may also be amplified and
   distributed to many types of downstream osmotic effects
   [[135]68,[136]69]. The cAMP as a signal factor participates in the
   regulation of osmotic pressure by mediating the increase of Na^+ uptake
   by cells in the internal environment [[137]70]. At the same time, the
   cAMP signaling pathway stimulates the production of arachidonic acid
   metabolites, which regulates the production of cortisol and glucagon,
   osmotic regulation and cellular fatty acid signaling in fish [[138]71].
   The “cAMP-responsive element modulator (CREM)” and “cAMP-specific 3′,
   5′-cyclic phosphodiesterase 4B (PDE4B)” were both down-regulated at 0
   psu and up-regulated at 45 psu. We suspected that the down-regulated
   expression of CREM and PDE4B at 0 psu inhibits Na^+ uptake and
   increases water transport, which maintains the balance of osmotic
   pressure. On the contrary, the up-regulation expression of CEMR and
   PDE4B at 45 psu enhances metabolic and osmolarity-regulated signals,
   providing energy and signaling mechanisms for maintaining osmolarity
   balance. Strangely enough, high affinity cAMP-specific and
   IBMX-insensitive 3′, 5′-cyclic phosphodiesterase 8A did decrease at 30
   psu, whether the specific mechanism of action related to the salinity
   tolerance of 30 psu is yet to be further explored. Furthermore,
   prolactin receptor (PRLR) was up-regulated at 0 psu for A. ommaturus.
   For freshwater adaptation, prolactin (PRL) is known to reduce Na^+
   efflux and water permeability, and is up-regulated under freshwater
   stress (0 psu) [[139]72], and PRL initiates its role by binding to a
   specific cell surface PRL receptor (PRLR) [[140]73]. In this way, PRLR
   participates in osmotic pressure regulation by affecting the binding of
   PRL hormones, and it indirectly regulates sodium ions to maintain the
   stability of osmotic pressure [[141]64]. It is predicted that all
   teleost species may have the PRL/PRLR system to adapt to sudden osmotic
   changes in the environment, and that PRLs are widely expressed in
   osmotically regulated organs [[142]64]. These signals have important
   functions in the signal cascade of osmotic regulation. It can be used
   to amplify or transmit osmotic sensory signals and then regulate the
   functional characteristics of various downstream osmotic effects.

5. Conclusions

   In this study, RNA-seq analysis was performed on A. ommaturus exposed
   to four salinity gradients (0 psu, 15 psu, 30 psu, and 45 psu), and the
   obtained high quality data could be the best guarantee of the
   reliability of experimental results. After the verification of this
   experiment, these changes are related to metabolic processes, ion
   transport, signal transduction and other characteristics. We also found
   that the similarity in the patterns of gene expression of gill ion
   transporters and channels in salinity adaptation between A. ommaturus
   and O. mossambicus. Compared with high salinity stress, A. ommaturus is
   better able to adapt to low salinity stress. This information will
   improve our understanding of salt tolerance of A. ommaturus at the
   molecular level and provide relevant transcriptome information for
   further study.

Supplementary Materials

   The following are available online at
   [143]https://www.mdpi.com/2073-4425/11/6/631/s1, Table S1: Primer
   sequences of the 8 target unigenes analyzed by qRT-PCR, Table S2:
   Summary of read statistics from RNA-sequencing of A. ommaturus, Figure
   S1: Top-hit species distribution, Figure S2: Proximal tubule
   bicarbonate reclamation, aldosterone-regulated sodium reabsorption,
   bile secretion, pancreatic secretion and mineral absorption.
   [144]Click here for additional data file.^ (2MB, pdf)

Author Contributions

   Conceptualization, N.S.; methodology, Z.S.; validation, Z.S. and F.L.;
   formal analysis, Z.S. and F.L.; resources, Y.Z. And N.S.;
   writing—original draft preparation, Z.S.; writing—review and editing,
   Z.S.; supervision, N.S.; project administration, N.S.; funding
   acquisition, N.S. All authors have read and agreed to the published
   version of the manuscript.

Funding

   This research was funded by National Key R&D Program of China
   (2018YFD0900905).

Conflicts of Interest

   The authors declare no conflict of interest.

References