RNA-seq and Analysis of Argyrosomus japonicus Under Different Salinities

Salinity variation affects the physiological processes of fish. This study analyzed the transcriptome of the gill tissue of Argyrosomus japonicus to determine the significantly differentially expressed genes (DEGs) of A. japonicus under salinity changes. Transcriptome analysis of nine samples yielded 55.873 Gb of clean data, 64,912 transcripts, and 29,567 unigenes, and 83.62% of the transcripts and 81.89% of the unigenes were annotated. Compared with the control group, the high- and low-salt groups showed 1,731 and 695 DEGs, respectively. Gene Ontology enrichment analysis revealed that the DEGs were significantly enriched in transportation, metabolism, and stress response. Kyoto Encyclopedia of Genes and Genomes pathway enrichment revealed that the DEGs were significantly enriched in some signaling pathways. Several key genes (KRT1, KRT2, ATP1A, LDH, PFN, ACTB_G1, TUBB, GZMB, MHC2, CCL19, EPX, ANXA5, ACBP, EHF, BHMT, COL1A, and RHOA) were related to salinity adaptation. When environmental salinity fluctuated, genes related to stress, immunity, ion transport, and metabolism became more sensitive. These results suggest that the adaptation of A. japonicus under salinity changes is a complex process that involves multiple genes acting together.


INTRODUCTION
Argyrosomus japonicus, a warm-temperature and near-bottom marine teleost fish of the order Perciformes, family Sciaenidae, and genus Argyrosomus (Chen and Zhang, 2015). The species resides in the nearshore sand and mud seabed area in the Yellow Sea, East China Sea, South China Sea, Southern Japan, Korean Peninsula, Indian Ocean, and Pacific Ocean (Chen and Zhang, 2015). It is an economically important marine fish in China, popular for its high nutritional and economic value (Han et al., 2021). A. japonicus is suitable for aquaculture because of its high price, marketability, high fecundity, fast growth, non-territorial or cannibalistic nature, and saline resilience (Fitzgibbon et al., 2007). Thus, this species is farmed in sea cages, coastal earthen ponds, and recirculating aquaculture systems (PIRSA, 2001).
The marine environment is changing rapidly around the world because of global warming (Philippart et al., 2011). Therefore, marine salinity is affected and changed. Salinity, as a key abiotic environmental factor, could affect the survival, growth and development, physiological activities, and immune metabolism of marine fish (Boeuf and Payan, 2001;Wang and Zhu, 2002;Kultz, 2015). Studies have shown some possible differences in the adaptability of different species to salinity changes, and similar species may have different responses to different degrees of salinity changes (Romano and Zeng, 2006;Ye et al., 2009). Several differentially expressed genes (DEGs) and pathways involved in salinity changes have been identified in many species, such as Oratosquilla oratoria (Lou et al., 2019), Rachycentron canadum , Crassostrea gigas (Zhao et al., 2012), Eriocheir sinensis , Exopalaemon carinicauda (Shen et al., 2020), and Nibea japonica (Meng et al., 2020). At present, more macroscopical studies on A. japonicus are available than those on the adaptive mechanism and the expression of functional genes and their regulatory pathways in response to salinity fluctuation.
Through the expression of different genes under salinity changes, some genes of salinity adaptability in A. japonicus could be found, and the salinity adaption mechanisms of A. japonicus could be understood. Transcriptome analysis provides an effective method to identify genes involved in salinity changes. Key genes and pathways that may be targets for selection can be identified by analyzing the transcriptional changes induced by salinity changes (Smith et al., 2013). In terms of aquatic fauna, transcriptome analysis has been successful in Ietalurus Punetaus , Oncorhynchus mykiss (Mohamed et al., 2012), and Penaeus orientalis . RNA-seq could also be used to understand the growth and development, adaptive evolution, immune response and stress response, and other biological processes of fish (Qian et al., 2014).
Studies on this species at the transcriptome level to identify the genes responsible for salinity regulation, which affects the understanding of the fundamental mechanism underlying adaptation to fluctuations in salinity, are lacking. In the present study, RNA-seq was performed on the gill tissues of A. japonicus under salinity gradient changes. Functional gene annotation, enrichment analysis, and differential expression were used to analyze the functional genes and biological processes involved in response to salinity changes. The results provided data for studying the molecular mechanism of salinity adaptation in marine fish. Furthermore, the RNA-seq data act as an important resource not only for the identification of novel genes but also for further investigations regarding A. japonicus biology.

Materials and Experimental Design
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.

RNA Extraction, Library Construction, and Illumina Sequencing
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 firststrand 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.

Transcriptome Quality Control, Assembly, and Functional Annotation
Raw RNA-seq data were filtered by removing reads with sequencing adapters, unknown nucleotides (N ratio > 10%), lowquality 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 1 alignments, with an E-value < 0.00001 (Shen et al., 2020).

Gene Differential Expression Analysis
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) 1 http://www.ncbi.nlm.nih.gov/BLAST/ FIGURE 1 | Argyrosomus japonicus.  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 2 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 Analysis
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 2 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.

Sequencing Data, Transcriptome Assembly, and Gene Annotation
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).

Gene Differential Expression Analysis
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 upregulated 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.

Co-expressed Differentially Expressed Genes
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 hyposalinity 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.

Gene Ontology Function Annotation and Enrichment Analysis of Differentially Expressed Genes
The  Figure 6C).

Validation of Transcriptomic Data by Quantitative Real-Time PCR
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 highthroughput 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.

DISCUSSION
In recent years, RNA-seq has been increasingly used in the study of transcriptome of many model and non-model plant and animal species, and remarkable results have been achieved (Hu Y. C. et al., 2015). With the continuous development of high-throughput RNA-seq, this technology has been widely used in the study of various aquatic animals. In the study of fish transcriptome, the main work is to screen DEGs and discover candidate genes. Transcriptomic analysis could probe gene expression and function at a genome-wide scale, which is likely to be useful for elucidating and analyzing the mechanisms underlying salinity stress adaptation in A. japonicus. In this study, the genes involved in osmoregulation environmental changes usually induce physiological responses in organisms, and the main response of organisms to external stress is gene expression (Evans and Hofmann, 2012;Wolf, 2013;. RNA-seq is a suitable method to identify the underlying molecular mechanisms associated with the regulation of osmotic pressure and immunity, and RNA-seq analysis related to osmotic pressure has been achieved in some aquatic vertebrates and invertebrates, such as Oreochromis mossambicus (Lam et al., 2014), Acipenser baerii (Guo et al., 2018), Oryzias latipes (Marty et al., 2014), Litopenaeus vannamei (Hu D. G. et al., 2015), E. sinensis , Ostrea gigas tnunb (Zhao et al., 2012), and Mytilus edulis (Lockwood and Somero, 2011). In the present study, the RNA-seq of A. japonicus showed that the differential gene expression levels in the transcriptome of A. japonicus under hyper-salinity were higher than those under hypo-salinity, indicating that more genes were involved in regulating the adaptation mechanism of A. japonicus exposed to hyper-salinity. The varied expression of genes indicated that many transcripts are differentially regulated in response to salinity acclimation. In principle, salinity changes may influence immunity, osmotic pressure, metabolism, and other physiological processes by changing the cellular physiological pathways. In fact, these cellular physiological pathways are regulated by some functional genes. Therefore, a close attention was paid to the functional genes associated with the osmotic pressure, immunity, and metabolism of A. japonicus exposed to salinity changes (Table 7 and Figure 9).

Genes Involved in Osmoregulation
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 TABLE 7 | Enriched DEGs may be related to osmotic adjustment, immunity, and energy metabolism.
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 upregulated 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.

Genes Involved in Immune Response
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.

Genes Involved in Energy Metabolism
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 hypersalinity, 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 longterm seawater (30 psu) acclimation (Rocha et al., 2005(Rocha et al., , 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 wholebody 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.

CONCLUSION
Genes related to stress, immunity, ion transport, and metabolism in the gill tissue of A. japonicus were sensitive to salinity changes. When the salinity of the external environment changes, the osmotic pressure could be regulated by adjusting the transport of substances and ions to maintain the balance of the internal environment. Salt resistance is not determined by a single gene; it requires many genes to work together. In this experiment, the transcriptome data of A. japonicus were obtained, and the sequence assembly, functional annotation, and various pathways of the expressed genes were analyzed. Genes related to salinity changes were obtained, which enhanced the transcriptome database of A. japonicus. The results could provide a scientific basis for the further study on the physiological mechanism of salinity changes and the in-depth study of key DEGs in A. japonicus. However, the specific molecular adaptation mechanism of A. japonicus under external salinity changes still needs to be further studied.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NBI SRA BioProject, accession no: PRJNA778869.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Care and Use Committee of Zhejiang Ocean University.

AUTHOR CONTRIBUTIONS
ZH and TG conceptualized the study and conducted the analyses. ZL and ZH analyzed the data and wrote and revised the manuscript. All authors contributed to the article and approved the submitted version.