Argyrosomus japonicus (Figure 1) was purchased from an aquatic product market in Zhoushan City, Zhejiang Province. This species was morphologically identified. Forty fish were temporarily kept in artificial seawater with a salinity of 25 psu. During the period, they were continuously oxygenated with oxygen pump, and the temperature was maintained at approximately 25°C. After 48 h, 30 randomly selected fish were exposed to different salinities for 24 h and stocked in three normal tanks (1.0 m × 0.7 m × 0.4 m) with 10 fish each. One group was exposed to the optimal salinity of 25 psu as the control group, and another two groups were exposed to the salinity of 35 and 15 psu as the high-salt group and the low-salt group, respectively. At the end of the experiment, three fish from each group were anesthetized, and samples were taken from the gill. The tissue was selected because it is considered a major osmoregulatory organ in most teleost fish (Higashimoto et al., 2001). These specimens were handled using RNA preservation solution and stored at a –80°C ultra-low temperature refrigerator for later use.

Total RNA was extracted from the gill tissues of the three groups by using a standard TRIzol Reagent Kit (Invitrogen, Carlsbad, CA, United States) in accordance with the manufacturer’s instructions. The integrity of RNA and the existence of DNA contamination were analyzed using 1% agarose gel electrophoresis. The absorbance values at the wavelengths of 230, 260, and 280 nm were measured using a NanoPhotometer spectrophotometer (Invitrogen, Carlsbad, CA, United States), and the ratios of OD260/280 and OD260/230 were calculated to detect RNA purity. RNA integrity was accurately measured using an Agilent 2100 bioanalyzer (Invitrogen, Carlsbad, CA, United States).

Illumina’s TruSeq RNA Library Preparation Kit was used following the manufacturer’s protocol for library construction. First, the mRNA with poly A tail was enriched using magnetic beads containing oligo (dT). Subsequently, the mRNA was randomly fragmented by adding NEB Fragmentation Buffer. The fragment mRNA was used as a template, and randomized oligonucleotides were used as primers to synthesize the first-strand cDNA. Then, the RNA chain was degraded by RNaseH. The second-strand cDNA was synthesized in the polymerase I system by using dNTPs as the raw material. The double-stranded cDNA was purified, and it underwent end repair. The poly A tail was added, and the sequencing joints were connected. The library was sequenced on Illumina HiSeq 2000, and the paired end was 150 bp.

Raw RNA-seq data were filtered by removing reads with sequencing adapters, unknown nucleotides (N ratio > 10%), low-quality reads (quality scores < 20), and non-AGCT bases in the 5′ end to ensure the accuracy of the subsequent bioinformatic analysis. The adapters and small fragments with a length less than 25 bp were removed after quality pruning. Thus, high-quality clean data could be obtained.

Comparison reference genome analysis was performed on Hisat2 (Kim et al., 2015; Johns Hopkins University), and the clean data obtained after quality control were compared with the reference genome of A. japonicus (GCA_015710095.1). Then, the mapped reads for transcript assembly and expression volume calculation were obtained. On the basis of the selected reference genome sequence, the software StringTie (Pertea et al., 2015) was used to assemble the mapped reads, compare them with the original genome annotation information, obtain the original unannotated transcription regions, and determine new transcripts and new genes of the species to supplement and improve the original genome annotation information. Sequences that coded for peptide chains that were too short (less than 50 amino acid residues) or contained only a single exon were filtered out. All unigenes were used to analyze the gene functions on the basis of NR (Non-Redundant Protein Sequence Database), Swiss-Prot (A manually annotated and reviewed protein sequence database), Pfam (Protein family database), STRING (An online search databases of known protein interaction relationships), GO (Gene Ontology), and KEGG (Kyoto Encyclopedia of Genes and Genomes) by using the BLAST (Scott and Madden, 2004) program¹ alignments, with an E-value < 0.00001 (Shen et al., 2020).

The expression levels of all unigenes were normalized, and RSEM and Bowtie2 were used to determine the FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) values (Mortazavi et al., 2008). The DEGs were identified across samples or groups by using DESeq2 (Love et al., 2014). FDR < 0.05 and | log₂FoldChange| ≥ 1.00 were considered as the filtering thresholds for significant DEGs. The selected DEGs were analyzed for the enrichment of GO and KEGG pathways, and P < 0.05 indicated that the DEGs were significantly enriched in the GO and KEGG pathways. P < 0.05 is generally considered as significant enrichment, and the smaller the P value is, the greater the significance.

Quantitative real-time PCR (qRT-PCR) was performed to verify the transcriptome data and the Illumina sequencing results. Five genes with a high level of significance were randomly selected in accordance with the sequencing data, and GAPDH (forward, 5′-AAGACCAACCCAGAGCAAAATG-3′; reverse, 5′-TCACCTTGAAGCGACCA-3′) was selected as the reference gene for qRT-PCR analysis. Gene-specific primers were designed using Primer Premier 6.0 (Premier Biosoft Inc., CA, United States), as shown in Table 1. Standard curves were constructed to determine the dilution times of cDNA samples, and they served as calibrators. qRT-PCR analysis was performed in accordance with the manufacturer’s instructions for ACEQ qPCR SYBR Green Master Mix (without ROX). A 10 μL reaction system was amplified by Bio-Rad CFX96 real-time PCR system (Applied Biosystems, Bio-Rad, United States), including 2 μL of cDNA, 5 μL of SYBR Green Mix (2×), 0.5 μL of forward primers, 0.5 μL of reverse primers, and 9.5 μL of ddH₂O. The PCR reaction conditions were as follows: 95°C for 5 min, 95°C for 10 s, and 60°C for 30 s, followed by 40 cycles. Dissolution curves were analyzed after the reaction was completed. Three parallel experiments were conducted for each cDNA template to strengthen the accuracy of the results. The relative expression of these genes was calculated using the 2^–ΔΔCt method as follows: ΔCt = Ct target unigene – Ct reference gene; ΔΔCt = ΔCt treatment – ΔCt control (Livak and Schmittgen, 2002). SPSS 19.0 statistical software (International Business Machines Inc., United States) was used for statistical analysis.

After filtering, approximately 55.873 Gb of clean data were generated for subsequent analysis. The percentage of Q20 base was above 97.50%, the percentage of Q30 base was above 93.70%, and the GC contents were all between 49.41 and 50.00%; the details are listed in Table 2. This finding indicated that the quality of RNA-seq of the gill tissue of A. japonicus improved. The clean data of each sample were aligned with the specified reference genome sequence, and the matching rates were above 92%, which indicated that the selected reference genome assembly could meet the requirements of information analysis.

A total of 64,368 transcripts were spliced into 29,567 unigenes after further redundancy removal (Table 3), and 40,802 new transcripts and 6,001 new unigenes were found. The target genes and transcripts were annotated, and 83.62% of the transcripts and 81.89% of the unigenes were annotated in the six large databases of NR, Swiss-Prot, GO, Pfam, KEGG, and STRING (Table 4).

The value of FPKM was taken as the gene expression amount in each sample. The expression levels of each sample are shown in a box diagram (Figure 2). The distribution of FPKM under different salinity conditions was not significantly different at the overall level. The correlation of gene expression levels between samples could be used to test the reliability of the experiment and whether the sample selection is reasonable. The sample correlation heatmap of gene expression (Figure 3) shows that the similarity of expression patterns among the samples was high, indicating high experimental reliability and reasonable sample selection.

The transcripts of the two experimental groups and the control group were analyzed. Figure 4A shows that the low-salt group had 695 DEGs, among which 305 were significantly up-regulated and 390 were significantly down-regulated. Figure 4B shows that the high-salt group had 1,731 DEGs, among which 847 were significantly up-regulated and 884 were significantly down-regulated. The number of significantly up-regulated and down-regulated DEGs in the high-salt group was higher than that in the low-salt group.

A Venn diagram was prepared to quantify the number of differentially expressed unigenes, and results showed that 22 unigenes were differentially expressed in the hyper- and hypo-salinity environments (Figure 5). As shown in Table 5, eight (EDC3; CYLD, USLP2; NCAM, CD56; PLA2G, SPLA2; KCNJ16; SLC9A3, NHE3; KRT1; E4.1.1.17, ODC1, speC, and speF) out of the 22 co-expressed genes were annotated to some specific KEGG pathways, which mainly involved RNA degradation, some immune-related signaling pathways, and substance metabolism and absorption. The result suggest that physiological processes, such as immunity response, signal transduction, substance metabolism, and osmotic pressure regulation, are part of the adaptation mechanism of A. japonicus to salinity fluctuations.

The main biological functions of DEGs could be determined by GO enrichment analysis. Under hyper-salinity, the GO enrichment results showed that the up-regulated DEGs were significantly enriched in transporter activity (GO: 0005215), transmembrane transporter activity (GO: 0022857), ion transmembrane transporter activity (GO: 0015075), chemical homeostasis (GO: 0048878), ion transport (GO: 0006811), transmembrane transport (GO: 0055085), anion transmembrane transporter activity (GO: 0008509), membrane part (GO: 0044425), intrinsic component of membrane (GO: 0031224), integral component of membrane (GO: 0016021), and multicellular organismal process (GO: 0032501) GO terms (Figure 6A). The down-regulated DEGs were significantly enriched in the DNA metabolic process (GO: 0006259), DNA repair (GO: 0006281), cell response to DNA damage stimulus (GO: 0006974), cell cycle process (GO: 0022402), and chromosomal part (GO: 0044427) GO terms (Figure 6C).

Under hypo-salinity, the up-regulated DEGs were significantly enriched in ion transport (GO: 0006811), inorganic molecular entity transmembrane transporter activity (GO: 0015318), transmembrane transporter activity (GO: 0022857), defense response to bacterium (GO: 0042742), oxoacid metabolic process (GO: 0043436), organic acid metabolic process (GO: 0006082), and positive regulation of secretion by cell (GO: 1903532) GO terms (Figure 6B); the down-regulated DEGs were significantly enriched in striated muscle contraction (GO: 0006941), extracellular matrix (GO: 0031012), and contractile fiber part (GO: 0044449) GO terms (Figure 6D).

Kyoto Encyclopedia of Genes and Genomes analysis showed that 14,724 unigenes were annotated and DEGs were enriched in 322 pathways. A large number of DEGs were involved in the regulation of the metabolic pathways (ko01100), pathways in cancer (ko05200), biosynthesis of secondary metabolites (ko01110), PI3K-Akt signaling pathway (ko04151), endocytosis (ko04144), MAPK signaling pathway (ko04010), HTLV-I infection (ko05166), neuroactive ligand-receptor interaction (ko04080), Ras signaling pathway (ko04014), proteoglycans in cancer (ko05205), Rap1 signaling pathways (ko04015), Focal adhesion (ko04510), cAMP signal pathway (ko04024), and regulation of actin cytoskeleton and axon guidance (ko04810), as shown in Table 6.

In accordance with the KEGG pathway enrichment analysis, the KEGG pathways involved in KEGG were divided into six categories: metabolism, cellular process, environmental information processing, genetic information processing, organismal systems, and human disease. Under hyper-salinity, the DEGs were mainly involved in circadian rhythm (ko04710), cGMP-PKG signaling pathway (ko04022), insulin resistance (ko04931), cAMP signaling pathway (ko04024), calcium signaling pathway (ko04020), and metabolic pathways (ko01100). When A. japonicus encountered hypo-salinity, the DEGs were mainly involved in the biosynthesis of amino acids (ko01230), proximal tubule bicarbonate reclamation (ko04964), IL-17 signaling pathway (ko04657), focal adhesion (ko04510), protein digestion and absorption (ko04974), rheumatoid arthritis (ko05323), ECM-receptor interaction (ko04512), and proteoglycans in cancer (ko05205), as shown in Figure 7.

Five target unigenes were selected to evaluate the transcriptome data of each experimental pair on the basis of qRT-PCR analysis. For these candidate unigenes, the variation trend in expression was identical between the qRT-PCR and transcriptome data, although the values derived from both analytical methods did not perfectly match (Figure 8). Consequently, the results indicated that the gene expression analysis based on high-throughput sequencing data was credible. The difference between the transcriptome and qRT-PCR results may be associated with the sequencing and mapping processes of transcriptome reads.

Salinity changes could directly affect the gene expression of osmotic pressure regulation. In this study, the pathways related to osmoregulation were mainly enriched when adapted to hyper-salinity. These pathways included proximal tubule bicarbonate reclamation (ko04964), aldosterone-regulated sodium reabsorption (ko04960), and mineral absorption (ko04978). ATP1A, ATP1B, GLS, NHE3, SNAT3, SLC4A4, SLC9A3R2, KCNJ1, CLCN2, and TRPV6 were significantly up-regulated, showing that they play a regulatory function in osmoregulation. When adapted to hypo-salinity, the proximal tubule bicarbonate reclamation (ko04964) related to osmoregulation was mainly enriched. SNAT3, NHE3, and AQP1 were significantly up-regulated, indicating that they participate in osmotic pressure regulation when salinity is reduced.

Sodium/potassium-transporting ATPase regulates K⁺ transmembrane input and Na⁺ transmembrane output to maintain osmotic pressure balance (Bloor et al., 2003). Na⁺/H⁺ exchanger, a transmembrane protein that exists in all vertebrate cells, is involved in the regulation of intracellular pH value, the control of cell volume, and ion transport (Pavel and Larry, 1998). Alexander et al. (2013) confirmed that NHE3 is involved in the reabsorption of Na⁺, bicarbonate, and water in the proximal tubules. Chloride channel participates in various physiological processes, including cell volume and membrane potential stabilization, ransepithelial transport, and signal transduction (Liu et al., 2015). Solute carriers, a family of membrane transport proteins, play important roles in physiological activities, such as material transportation, energy transfer, and signal transduction between cells (Wang et al., 2017). Aquaporins constitute a small hydrophobic family that allows osmotic driving water and small solutes to pass through biofilms, and they are found in all living organisms (Calvanese et al., 2013). By adjusting the plasma concentrations of Na⁺, K⁺, and Cl^– inside and outside the cell and the reabsorption of Na⁺, bicarbonate and water, A. japonicus adjusts its own osmotic pressure to maintain its balance with the fluctuation of salinity.

Analysis showed that the calcium signaling pathway was also enriched, and certain genes related to the calcium signaling pathway were identified. In particular, ADRB2, HRH1, AGTR1, EDNRA, OXTR, TBXA2R, CACNA1I, and VDAC2 were up-regulated under hyper-salinity. The calcium signaling pathway is the most researched signaling pathway at present. Fiol et al. (2006) proposed that the calcium signal transduction in fish gills is an important signal transduction mechanism for fish gills to respond to osmotic stress. Marshall et al. (2000) confirmed the importance of the calcium signaling pathway in fish corresponding to salinity changes. However, the specific signal transduction mechanism in the regulation of osmotic pressure needs to be further studied.

In the present study, the genes related to ornithine transcarbamylase, argininosuccinate synthetase, and carbamoyl phosphate synthetase-I were up-regulated at hypo-salinity, and they regulated osmotic pressure by participating in urea production.

Salinity changes disrupt the immune system, causing weakness in the organism defense of bacteria and diseases. Immune responses are important in adaptation to salinity changes. The present research showed that DEGs were enriched in some immune pathways under hypo-salinity, including IL-17 signaling pathway (ko04657) and cytokine–cytokine receptor interaction (ko04060). The related IL1B, IL17C, MUC5AC, IL22RA2, and CCL19 were up-regulated.

Fish are ectotherms, and their innate immune response mechanism plays a major role in immune response (Bai et al., 2017). Some immune-related cytokines were associated with salinity reduction. Cytokines are only secreted in large quantities when the body enters the stage of immune response (Secombes et al., 1996). The cytokines identified in fish are interleukins (ILs), chemokines, interferon, and transforming growth factor-β (Bai et al., 2017). IL plays an irreplaceable role in fish innate immunity and adaptive immunity, and the pro-inflammatory activities of IL-1 cytokines are crucial during the activation of the innate immune system (Iwasaki and Medzhitov, 2015; Bai et al., 2017). Chemokines could activate chemotactic leukocytes in inflammation (Baoprasertkul et al., 2005). C-C chemokines are an important part of the innate immune system, and they recruit white blood cells and enhance innate immune response (Hu and Zhang, 2015).

Given their high expression, some cytokines and chemokines at hypo-salinity activate the innate immune system of A. japonicus to suppress inflammation and disease. All immune-related genes could work together to trigger complex defense mechanisms.

The mechanisms involved in environmental stress are energy consuming, and they ultimately influence the future life characteristics of the response of organisms (Verslycke and Janssen, 2002). When salinity changes, osmotic pressure regulation is an energy-consuming process. Therefore, an essential prerequisite for successful salinity adaptation is the activation of energy metabolism. In the present study, some pathways related to energy metabolism were significantly enriched under hyper-salinity.

In multicellular organisms, ATP, the major energy source in the body, is predominantly supplied by a series of metabolic pathways, including glycolysis, the citric acid cycle, and the electron transport chain (Nelson and Cox, 2008). Glycolysis occurs with variations in nearly all organisms, aerobic and anaerobic. The wide occurrence of glycolysis indicates that it is an important energy metabolic pathway (Romano and Conway, 1996). In the present study, glycolysis/gluconeogenesis (ko00010) related to energy metabolism was enriched. The present study showed that the ACSS, PK, LDH, and ALDO related to glycolysis were up-regulated. Lactate dehydrogenase (LDH) is an important enzyme for glycolysis in anaerobic environment. Pyruvate kinase (PK) regulates ATP, ADP, and glycolysis intermediate products in cells, and it plays a key role in the control of glycolysis rate (Wang et al., 2005). Acetyl-coenzyme A synthetase provides the cell with two carbon metabolites used in many anabolic and energy generation processes (Starai and Escalante-Semerena, 2004). ACSS catalyzes the ligation of acetate with CoA to produce acetyl-CoA, which plays important roles in fatty acid and cholesterol synthesis and tricarboxylic acid cycle (Fujino et al., 2001). Thus, to maintain the osmotic and ionic homeostasis under hyper-salinity, the high expression levels of ACSS ensure that more acetyl-CoA from energy substance metabolism can enter the Krebs cycle and promote the reaction continuously. In addition, the high expression levels of LDH and PK promote glycolysis and provide more energy support.

In the present study, the genes with high expression may be related to lipid metabolism. For example, adipose triglyceride lipase is a major rate-limiting enzyme that regulates lipid storage and release in adipocytes (Miyoshi et al., 2008). Previous works discussed lipid metabolism in terms of fish osmoregulation because lipids are an important energy source in fish. Previous studies demonstrated the importance of lipids for meeting the requirement of metabolism in Centropomus parallelus after long-term seawater (30 psu) acclimation (Rocha et al., 2005, 2007). The DEGs were also enriched in the AMP-activated protein kinase (AMPK) signaling pathway (ko04152). AMPK is allosterically activated by AMP, and AMPK is necessary to regulate whole-body energy balance (Hardie and Sakamoto, 2006; Towler and Hardie, 2007). The present study showed that A. japonicus needs energy to regulate a series of physiological changes caused by hyper-salinity.