1 INTRODUCTION

The kuruma shrimp Marsupenaeus japonicus is one of the most commercially important shrimp species worldwide, distributed in Japan, the South China Sea, and the Indo-West Pacific (Hamasaki and Kitada, 2008; Tsoi et al., 2014). The adaptability of M. japonicus to salinity is relatively poor, the suitable salinity range is 15-30, the survival salinity is 7-35, and it is difficult to survive below 7 (Weng et al., 2012). Due to its advantages of rapid reproduction, tolerance to low temperature and long-distance transport capacity, this species has become one of the most important shrimp species in China's coastal waters (Paterson, 1993; Duan et al., 2016). In 2016, the production of M. japonicusin China was 55 885 t (http://www.fao.org/home/en/). In recent decades, there have been many studies on M. japonicus, and mostly concerned about the life activities (Chen and Lai, 1993; Hewitt and Duncan, 2001), cultivation mode (Cheng and Chen, 2001; Li et al., 2014b), environmental stress (Coman et al., 2002; Lee and Chen, 2003), and molecular biology (Lin et al., 2008; Leu et al., 2011). Recently, production of M. japonicus has declined; in particular, production was approximately 2 000 t less in 2017 than in 2016 (http://www.fao.org/home/en/). This decrease may be related to the broodstock of M. japonicus collected from wild populations, which exhibit poor stress resistance (Liu et al., 2019). The another, it may be the environmental stress, which combined with a bottom environment rich in ammonia nitrogen and nitrite nitrogen in farms undergoing a series of environmental changes (Cheng et al., 2004, 2013), such as salinity (Setiarto et al., 2004), temperature (Lu et al., 2016), and dissolved oxygen (Chien and Shiau, 2005).

Salinity is one of the most important environmental factors affecting the survival of aquatic organisms (Heugens et al., 2001; Hop et al., 2002; Laing, 2002; Chong-Robles et al., 2014). With the rapid change of the global marine environment, ocean warming (Johnson et al., 2011), typhoon activity (Balaguru et al., 2016), tide movement (Kourafalou et al., 1996) and other natural and human factors can cause the decrease of the salinity of the natural sea water, which has a great impact on aquatic organisms. Low salinity stress has been reported in many shrimps, such as Litopenaeus vannamei (McGraw et al., 2002), Fenneropenaeus chinensis (Liu et al., 2006), Penaeus monodon (Ye et al., 2009), and M. japonicus (Via, 1986). Low salinity affected the growth and survival of prawn, and the survival rate decreased significantly with the decrease of salinity (Via, 1986; McGraw et al., 2002; Ye et al., 2009).In addition, low salinity affected osmotic regulation (McNamara et al., 2015; Fregoso-López et al., 2017) and immune regulation (Lin et al., 2012). For example, M. japonicus and L. vannamei are susceptible to white spot virus and Vibrio alginolyticus under low-salinity conditions (Yu et al., 2003; Wang and Chen, 2005; Li et al., 2010), which may be due to shrimp consuming a large amount of stored energy to regulate osmotic pressure in vivo, leading to an increase in metabolic rate and decreases in growth rate and stress tolerance upon salinity change (Ingram et al., 2002; Setiarto et al., 2004). To clarify the osmotic regulation mechanism used by shrimp under low-salinity conditions, molecular research has intensified, including studies using gene cloning (Tiu et al., 2007; Meng et al., 2011; Rajesh et al., 2012), suppression subtractive hybridization (Gao et al., 2012; Shekhar et al., 2013), transcriptome sequencing (Hu et al., 2015), proteomics (Fan et al., 2019), and digital gene expression techniques (Zhao et al., 2015). Among these methods, transcriptome sequencing has been used most widely. Unfortunately, there is still a lack of research on low salt molecules of M. japonicus, the specific mechanism of osmotic and immune regulation remains unknown.

Transcriptome sequencing can be used for accurate quantification of the expression of various RNA molecules, such as mRNAs, lncRNAs, and microRNAs, in different tissues, developmental stages, and disease states, to characterize complex molecular pathways and differential gene expression patterns (Brawand et al., 2011; He et al., 2012; Zhao et al., 2012; Park et al., 2013; Kornienko et al., 2016; Ludwig et al., 2016). To date, the transcriptomes of several prawn species, including Litopenaeus vannamei (Li et al., 2012), Fenneropenaeus chinensis (Li et al., 2013), Fenneropenaeus merguiensis (Powell et al., 2015), and Macrobrachium rosenbergii (Jung et al., 2011), have been sequenced. Those studies focused mainly on development, immunity, and osmotic regulation in shrimp. Previous studies have shown that the main osmotic regulator of shrimp is the gill (Shekhar et al., 2013; Hu et al., 2015). However, recent studies have found that the hepatopancreas is also involved in osmotic regulation (Chen et al., 2015). The hepatopancreas is the main organ responsible for not only the digestion, metabolism, absorption, and storage of nutrients in decapod crustaceans but also for energy metabolism under low-salinity stress (Wang et al., 2014b; Zhang et al., 2016). The purpose of this study is to evaluate transcript levels in the hepatopancreas of M. japonicus under low-salinity stress for different periods using transcriptome sequencing, identifying the pathways related to osmotic regulation, immune response, detoxification, and apoptosis, and revealing differentially expressed genes (DEGs) involved in low-salinity stress. The results not only enrich the information of shrimp cDNA database, but also provide important molecular evidence for breeding improved strains of M. japonicus in the future.

2 MATERIAL AND METHOD 2.1 Salinity challenge and sample preparation

A 120-day-old mixed family of M. japonicus with an average body length of 49.28±4.79 mm and average body weight of 1.39±0.38 g was obtained from Guolian Aquatic Limited Company, Guangdong Province, China. Before the experiment, the shrimp were acclimated in a cement pond in 28±0.2 ℃, salinity 29.8±0.2, and dissolved oxygen 6.0 mg/L for 7 days to reduce any stress response. In the preliminary experiment, the adaptability of M. japonicas to low salinity was studied. The results show that the survival rate of M. japonicus was 80% for 96 h in salinity 8. Thus, we determined that the salinity of low salt stress was 8 in this study. The shrimp were divided into the control group (CG) and low-salinity group (LS). Salinity was regulated by mixing seawater with disinfected tap water, and the low-salinity and control groups were exposed to salinities of 8 and 29.8±0.2, respectively. Three experimental replicates were performed per group. A 1 500-L plastic bucket was used for the experiment, and all other experimental conditions were consistent with the domestication conditions. No food was given during the experiment, which lasted 96 h, with salinity determined every 6 h to maintain a constant level. The hepatopancreas was collected from 9 shrimp of each group at the 6^th, 12^th, 24^th, 48^th, and 96^th h and stored in an enzymatic centrifuge tube containing 1-mL RNA stabilization solution. The samples were stored overnight at 4 ℃ and then held at 20 ℃ until RNA extraction.

2.2 RNA extraction and Illumina sequencing

Total RNA from the hepatopancreases was extracted using TRIzol reagent (Invitrogen) according to the manufacturer's instructions. RNA degradation and contamination were monitored by 1% agarose gel electrophoresis. RNA purity was checked using the NanoPhotometer^® spectrophotometer (Implen, CA, USA). The RNA concentration was measured using the Qubit^® RNA Assay Kit and Qubit^®2.0 fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) (Hu et al., 2015). Ten libraries were constructed, five from the low-salinity group at each time point (LS6, LS12, LS24, LS48, and LS96) and five from the control group at each time point (CG6, CG12, CG24, CG48, and CG96). These 10 libraries were generated using the NEBNext^® Ultra^TM RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer's recommendations, and index codes were added to identify the sequences expressed in each sample. The total RNA content of each sample was 3 μg. Finally, PCR products were purified (The Agencourt AMPure XP system, Beckman Coulter, Brea, CA, USA), and the library quality was assessed on the Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was performed on the cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, library preparation was carried out on the Illumina Hiseq platform, and 125–150 bp paired-end reads were generated (Fu and Zou, 2018).

2.3 Bioinformatics analysis of the transcriptome data 2.3.1 Sequencing data filtering and assembly

Raw data (raw reads) in the fastq format were first processed using in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter sequences or poly-N sequences and low-quality reads (The base number of the qphred≤20 accounts for more than 50% of the total reads) from the raw data. At the same time, the Q20, Q30, and GC content of the clean reads were calculated, respectively. All downstream analyses were based on high-quality clean reads (Li et al., 2020). Clean reads were assembled with Trinity (Grabherr et al., 2011). Briefly, Trinity is an efficient and stable transcriptome splicing software developed by Broad Institute and Hebrew University of Jerusalem for RNA-seq data. It combines three independent software modules (Inchworm, Chrysalis, and Butterfly) to process and splice a large number of RNA-seq. The main process was as follows: By software Inchworm, first, all the fq files of the read were read and converted into the fa format read files, and then merged the 3′-end reads and 5′-end reads to obtain both.fa. The reads were decomposed by kmers (short fragment of k-bp length), and the types of k-mers and the number of k-mers of each type were counted, and then sorted from high to low by the frequency of k-mers. The k-mer with the highest frequency, as a seed, was extended to 3′-end. The frequency of each k-mer after the extension was counted. The k-mer with the highest frequency was selected as the extension path. Finally, following the greedy algorithm, the k-mers were extended by the overlap relationship to form the contig sequence. Using software Chrysalis, all contigs with similar regions more than k–1-mers were clustered to form components. The de Bruijn graph was constructed according to different components, and meantime, the reads and components were compared and verified. With the software Butterfly, the transcripts were listed, and the graph was split into linear sequence. The error sequence was eliminated by the read and pairs relationship, and the de Bruijn graph of each component was simplified. The full-length transcripts of the variable splicing subtypes were output, and the transcripts corresponding to the paralogous genes were sorted out. Finally, the splicing result file is obtained and named as TRINITY. fasta.

2.3.2 Annotation of gene functions and differential expression analysis

Gene functions were annotated based on the following databases: National Center for Biotechnology Information non-redundant protein (NR) and nucleotide sequences (NT), Pfam (protein family database), Clusters of Orthologous Groups of proteins (KOG/COG), Swiss-Prot (a manually annotated and reviewed protein sequence database), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2007), KEGG Orthology database (KO) (Mao et al., 2005), and gene ontology (GO) databases (Young et al., 2010).

Clean reads were aligned to transcriptome sequences by bowtie2 (Langmead et al., 2009), and the alignment results of bowtie were counted by RSEM (http://deweylab.biostat.wisc.edu/rsem/) (Li and Dewey, 2011). The read counts of each sample aligned to each gene were obtained, and the FPKM conversion were carried out (Trapnell et al., 2010). Read counts for each sequenced library were adjusted using the edgeR v.3.0.8 package of R software and a normalization scaling factor, TMM was used to standardize the read count data, and the DEGs between LS groups and CG groups were analyzed with DEGseq R package. Differential genes were determined by combining q value (q value≤0.05) and fold change (|log₂FC(sample2/sample1)|≥1). P values were adjusted using q values. A q value of 1 was set as the threshold for defining significant differential expression (Wang et al., 2010).

2.4 Validation of mRNA-seq data using quantitative real-time PCR (qPCR)

Nine DEGs from five libraries were randomly selected, including taurine catabolism dioxygenase (TauD), sarcosine dehydrogenase mitochondrial isoform X1 (SarDH), C-type lectin (CLEC), glycine N-methyltransferase (GNMT), betaine-homocysteine S-methyltransferase (BHMT), methionine adenosyltransferase II alpha (MAT2α), serine/ threonine-protein phosphatase 2A catalytic subunit beta isoform (PPP2Cβ), taurine transporter (TAUT), and agmatinase mitochondrial-like (AGMAT) were confirmed by qPCR. Using TransStart Tip Green SuperMix (Beijing TransGen Biotech Co., Ltd.) and performed on the Light Cycler 480II (Roche) thermal cycler according to the manufacturer's instructions. Elongation factor 1α (EF1α) was used as the reference gene, and each gene amplification was performed in three wells. Gene-specific primers were designed using Premier 5 (Table 1). Amplification was performed in a 384-well plate in a 10-μL reaction volume containing 1-μL cDNA, 0.2 μL each of the gene-specific forward and reverse primers, 5-μL TransStart Tip Green qPCR SuperMix and 3.6-μL RNase-free water. The thermal profile for the PCR was 94 ℃ for 30 s followed by 45 cycles of 94 ℃ for 5 s, 60 ℃ for 15 s and 72 ℃ for 10 s. The results are presented as the change in relative expression normalized to the reference gene (EF1α) using the 2^-ΔΔCt method. The data was mean±SD, and the difference was statistically significant (P < 0.05) by one-way ANOVA in SPSS19.0.

3 RESULT 3.1 Sequencing and de novo assembly of the transcriptome

After transcriptome sequencing, 92 184 260– 108 448 152 raw reads were generated from the 10 libraries prepared from the control and low-salinity groups. After removing the reads containing adapter sequences or poly-N sequences and low-quality reads from the raw reads, 88 890 960–1 051 300 444 clean reads were obtained, the data volume is 6.67G–7.88G, the GC contents of the raw data were 50.29%–52.48%, the Q20 percentages were 97.62%–98.11%, and the Q30 percentages were 93.34%–94.53% (Table 2). Using the Trinity to assemble the obtained clean reads, and after removing the redundancy, 48 807 unigenes were obtained from 52 826 925 nucleotides. N50 and N90 were 1 746 bp and 436 bp, respectively, with an average length of 1 082 bp (Table 3). According to the length distribution statistics of the assembled unigenes, the minimum length is 301 bp; there are 19 581 unigenes distributed between 300– 500 bp, accounting for 40.12% of the total, with the largest number; there are only 6 383 unigenes more than 2 000 bp, accounting for 13.08% of the total (Table 4; Fig. 1).

We used BUSCO to evaluate the splicing quality of unigenes, and compared with the conservative genes to evaluate the accuracy and integrity of the splicing results. The results showed that 978 BUSCO were completely covered, and 918 (93.9%), 15 (1.5%), 20 (2.0%), and 25 (2.6%) unigenes were completely matched, including complete and single copy, complete and duplicated, partial fragment, and missing, respectively (Table 5; Fig. 2). These data greatly enrich the genetic resources available for M. japonicus, which may facilitate further molecular studies of M. japonicus under low-salinity conditions for various periods.

3.2 functional annotation of genes

A total of 48 807 unigenes were annotated by NR, NT, KO, Swiss-Prot, Pfam, GO, and KOG. The number of unigenes annotated to at least one database in seven databases is 36 707, accounting for 75.20% of the total, while the number of unigenes annotated to all seven databases is 1 466, only accounting for 3.00% of the total (Table 6). Among them, in the light of the analysis of species sources, 31.6%, 14.4%, and 6.6% of the transcript data compared to the NR database were annotated to the top 3 species, Hyalella azteca, M. japonicus, and Litopenaeus vannamei (Fig. 3).

According to the GO annotation of unigenes, 26 763 unigenes obtained GO annotation, accounting for 54.83% of the total. The annotation can be divided into three categories: biological process (BP), molecular function (MF), and cell component (CC). In the biological process, the most involved are cell process (56.80%), metabolism process (52.31%), and single cell tissue process (49.19%). In the cell component, cell part, cell, and membrane account for 26.49%, 26.49%, and 20.20% of the total, respectively. In molecular function, it is mainly related to binding and catalytic activity, accounting for 47.56% and 46.31% of the total, respectively (Fig. 4).

A total of 9 919 (20.32%) unigenes were successfully annotated in the KOG database and classified according to 26 KOG, among which the most unigenes annotated to R (general function prediction) accounted for 14.30% of the total, followed by O (post translation modification, protein transformation, and molecular chaperone), T (signal transduction mechanism), E (amino acid transport and metabolism), accounting for 9.87%, 9.17%, and 8.56% of the total, respectively. However, X (unnamed protein) accounts for only 0.02% of the total, with the least annotated unigenes (Fig. 5).

After KO annotation of 3 536 (7.24%) unigenes, according to the KEGG metabolic pathway involved by unigenes, there were 485 (13.72%), 427 (12.08%), 553 (15.64%), 2 006 (56.73%), and 795 (22.48%) unigenes involved in the A, B, C, D, and E (cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems) of KEGG metabolic pathway, respectively Middle branch. Among all pathways, the most abundant unigenes were 395 (11.17%) in carbohydrate metabolism. Secondly, there are 351 and 337 unigenes in amino acid metabolism and signal transduction, respectively, accounting for 9.93% and 9.53% of the total (Fig. 6).

3.3 Differential gene expression analysis

The annotated genes in the control and low-salinity groups were screened to obtain DEGs. A total of 811, 589, 1 095, 745, and 875 DEGs were identified in the LS6, LS12, LS24, LS48, and LS96 libraries compared with the CG6, CG12, CG24, CG48, and CG96 libraries, respectively (Fig. 7). More genes were upregulated than down-regulated among all of the DEGs in 6 h, with approximately double the number of upregulated than down-regulated genes. However, the number of down-regulated genes exceeded the number of up-regulated genes in 12 h, 48 h, and 96 h.

3.4 Functional annotation of DEGs

For DEG annotation, a total of 811, 589, 1 095, 745, and 875 genes in the LS6, LS12, LS24, LS48, and LS96 groups, respectively, were identified using the NR, NT, KO, Swiss-Prot, Pfam, and KOG databases (Table 7). The percentage is ranked as NR > Pfam > Swiss-Prot > KOG > NT > KO, of which the percentages of NR were all greater than 80%, while those of the KO database were 9.34%–11.89%.

3.5 GO and KEGG enrichment analysis of DEGs

The DEGs from the comparisons of LS6 vs. CG6, LS12 vs. CG12, LS24 vs. CG24, LS48 vs. CG48, and LS96 vs. CG96 were analyzed based on GO enrichment and classified according to the three major GO categories: biological processes, cellular components, and molecular functions. There are 1 822–2 388 GO terms types. DEGs are the most abundant in biological processes, accounting for 56.97%–60.29% of the total, followed by molecular function (Table 8). Among the numerous biological processes represented, metabolic processes were the best represented. Chitin metabolic process (GO: 0006030), glucosamine-containing compound metabolic process (GO: 1901071), and amino sugar metabolic process (GO: 0006040) were the most commonly enriched among the groups. Only the LS6 vs. CG6, LS24 vs. CG24, and LS48 vs. CG48 comparisons showed significant enrichment of cellular component-related DEGs, with most related to extracellular components and the nucleosome. In addition, genes related to chitin and iron binding were significantly enriched in all groups among molecular functions (Fig. 8).

DEGs were mapped to KEGG pathways to confirm the biological pathways activated or repressed during the response of M. japonicus to low salinity. A total of 151, 85, 149, 97, and 111 KEGG pathways were enriched among the DEGs identified in the LS6 vs. CG6, LS12 vs. CG12, LS24 vs. CG24, LS48 vs. CG48, and LS96 vs. CG96 comparisons, respectively, among which the first 20 significantly enriched pathways are shown (Fig. 9). The results show that metabolic pathway and immune pathway were significantly enriched, which may play important roles related to salinity stress include the Linoleic acid metabolism (ko00591), Arachidonic acid metabolism (ko00590), Sphingolipid metabolism (ko00600), AMPK signaling pathway (ko04152), Sphingolipid signaling pathway (ko04071), Lysosome (ko04142), and TGF-beta signaling pathway (ko04350).

3.6 Verification of the significant DEGs by qPCR

Nine DEGs were selected to verify the RNA-seq results in the comparisons of the experimental (LS6, LS12, LS24, LS48, and LS96) and control (CG6, CG12, CG24, CG48, and CG96) groups using qPCR. The qPCR results showed that the trends in DEG variations were consistent with the RNA-seq results, verifying the accuracy of RNA-seq (Fig. 10). According to qPCR, TauD, SarDH, CLEC, GNMT, BHMT, MAT2α, and PPP2Cβ showed up-regulation in the control group, whereas TAUT and AGMAT showed down-regulation in the control group.

4 DISCUSSION 4.1 Ion exchange-related processes under salt stress

Decapod crustaceans exhibit a wide range of osmoregulatory patterns and capabilities (McNamara and Faria, 2012). In the low ionic freshwater environment, decapods face four major osmoregulation challenges simultaneously: absorbing ions from food and the surrounding water, minimizing diffusive ion loss via semi-permeable membranes, restricting incoming water flux and regulating cell volume decreases to maintain cellular stability. These osmotic and ionic regulation processes are conducted mainly in the gills (Freire et al., 2008; Furriel et al., 2010; Rahi et al., 2017). Numerous enzymes and transporters are involved in ion transport in crustacean gills, including Na⁺/K⁺-ATPase, V (H⁺)-ATPase, carbonic anhydrase and Na⁺/K⁺/2Cl cotransporter (Weihrauch et al., 2001; Tsai and Lin, 2007; Furriel et al., 2010; Havird et al., 2013). Interestingly, we found three biological processes related to osmotic regulation in the LS12 vs. CG12 comparison of M. japonicus under low-salinity stress: energy-coupled proton transmembrane transport (GO: 0015988), ATP hydrolysis-coupled transmembrane transport (GO: 0090662) and ATP hydrolysis-coupled transport proton (GO: 0015991). These three processes were enriched among the DEGs identified at the 12^th h only, and almost all of these DEGs were down-regulated. This result shows that the hepatopancreas of M. japonicus was involved in osmotic regulation, via inhibition of the ion transport process, under low-salinity stress for 12 h. We found significant differences in two DEGs: V (H⁺)-ATPase and V-H-ATPase subunits. V (H⁺)-ATPases are a family of membrane-related ATP-dependent proton pump enzymes. When located in the plasma membrane, these proteins couple the energy released during ATP hydrolysis with active transport of protons from the cytoplasm to intracellular or extracellular spaces (Forgac, 2007). The V1 region of V (H⁺)- ATPases contains eight different subunits (A–H), with catalytic sites located at the interface between the A and B subunits participating in ATP hydrolysis (Jefferies et al., 2008). When L. vannamei was moved from seawater with a salinity of 20 to water with a salinity of 5, the expression of V-H-ATPase subunit A in the gills was increased at approximately 12 h after the change in salinity, whereas that of V-H-ATPase subunit A in the hepatopancreas was decreased (Wang et al., 2012). In this study, we found that the expression levels of V(H⁺)-ATPase and V-H-ATPase subunits were down-regulated in 12 h, similar to previous findings in L. vannamei. We speculate that the hepatopancreas of M. japonicus inhibits H⁺ transport from the intracellular to extracellular space as an adaptation to hypoosmotic stress. The specific reason for this response remains to be explored.

The ion transporter Na⁺/K⁺-ATPase alpha subunit was also differentially expressed in the LS24 vs. CG24 comparison. Ion absorption by crustaceans is achieved mainly via the Na⁺/K⁺-ATPase (Lucu and Towle, 2003). The Na⁺/K⁺-ATPase is located in the basolateral membrane of gill ion-transport epithelial cells. The establishment of an electrochemical gradient shows that the Na⁺/K⁺-ATPase is the main force driving ion transport (Towle and Kays, 1986). The Na⁺/K⁺-ATPase maintains a gradient of excess Na⁺ outside and excess K⁺ inside the cell, and this gradient is used as an energy source to create a membrane potential for regulating cytoplasmic composition and transepithelial transport (Skou and Esmann, 1992; Suhail, 2010). Expression of the Na⁺/ K⁺-ATPase alpha subunit gene in the hepatopancreas was increased under hypoosmotic stress, but the transcript level of this gene was down-regulated in 24 h, which was not consistent with previous research (Li et al., 2015). This pattern may result from the differences of species, and the specific reasons require further exploration.

4.2 Lipid metabolism and glucose metabolism related to osmotic regulation under salt stress

Under low-salinity stress, shrimp require additional saturated and polyunsaturated fatty acids to maintain their osmotic regulation and ion exchange, and this process involves high-energy consumption (Chen et al., 2014, 2019). A number of lipid metabolic pathways in the hepatopancreas of L. vannamei are affected by salinity stress, including fatty-acid biosynthesis, arachidonic acid metabolism, adipocytokine signaling, and glycerophospholipid metabolism (Chen et al., 2015). Fatty acid metabolism and arachidonic acid metabolism are also related to osmotic regulation in Eriocheir sinensis (Li et al., 2014a). In this study, enrichment of DEGs associated with linoleic acid metabolism (ko00591) and arachidonic acid metabolism (ko00590) was found in the comparisons of LS6 vs. CG6, LS24 vs. CG24, LS48 vs. CG48, and LS96 vs. CG96, and most of these DEGs were up-regulated. We speculate that lipid metabolism in the hepatopancreas (the main organ of fat production) and lacrimal secretion are increased in M. japonicus. Acceleration of synthetic renewal can provide enough energy for osmotic regulation in shrimp at low salinity. Previous studies have revealed two major strategies for osmotic regulation in crustaceans, the "limiting process" and "compensatory process", which are carried out mainly by the gills. The limiting process is a strategy to maintain hemolymph osmotic pressure and ions in branchial membranes by adjusting the permeability of boundary structures, which can effectively reduce ion diffusion and water inflow (Péqueux, 1995; Rainbow and Black, 2001). Our results suggest that the strategy of the limiting the process may play an important role in osmotic regulation in M. japonicus, which is similar to that in Portunus trituberculatus (Lv et al., 2013) and L. vannamei (Chen et al., 2019).

Energy metabolism plays an important role in survival and normal body functions, as well as in stress adaptation and tolerance (Sokolova et al., 2012). In this study, we found that glycometabolism was significantly affected by low-salinity stress in M. japonicus, with changes in fructose and mannose metabolism (ko00051), sphingolipid metabolism (ko00600), and other glycan degradation (ko00511) pathways and significant up-regulation of the glucosamine-containing compound metabolic process (GO: 1901071), amino sugar metabolic process (GO: 0006040), and aminoglycan metabolic process (GO: 0006022) terms in all groups. DEGs related to glycometabolism in the hepatopancreas of L. vannamei were more represented at a salinity of 3 than in the control in terms of gluconeogenesis, amino sugar and nucleotide sugar metabolism, glycosphingolipid biosynthesis, glycosaminoglycan degradation, and other glycan degradation, suggesting a requirement for extra energy by L. vannamei under chronic low-salinity stress, which is usually derived from the metabolism of nutrient substances (Xu et al., 2017). Therefore, we speculate that the active glycometabolism of M. japonicus under low-salinity stress is an adaptation to maintain osmotic regulation. Fructose-1, 6-bisphosphate aldolase is a key enzyme involved in glycolysis and gluconeogenesis. It is ubiquitously expressed in animals, plants and microorganisms (De-Simone et al., 2006; Lu et al., 2012; Dawson et al., 2013). In this study, fructose- 1, 6-bisphosphate aldolase was up-regulated at the 6th h and then returned to a normal level, reflecting an increase in glycometabolism during the early stage of low-salinity stress in M. japonicus. Previous studies have shown that supplementation with appropriate carbohydrates can improve the stress resistance of L. vannamei to low salinity (Wang et al., 2014c). Carbohydrate metabolism appears to play a major role in the energy supply for iono- and osmoregulation, and the liver is the main source of carbohydrate metabolites for osmoregulatory organs (Tseng and Hwang, 2008). Up-regulation of related genes indicated that carbohydrate metabolism was increased to provide more energy during the early stage of low-salinity stress in M. japonicus. We identified a DEG, alpha-amylase (alpha-1, 4-glucanglucanohydrolase), that was significantly associated with carbohydrate metabolism. Alpha-amylase belongs to the glycosyl hydrolase family 13, and it is a major carbohydrate hydrolase involved in the digestive functions of the liver and pancreas in shrimp (Kuriki and Imanaka, 1999; Glenn et al., 2005). In this study, alpha-amylase expression returned to normal after 6 h of up-regulation and then was downregulated after 24 h. Therefore, we speculate that this gene is related to low-salinity stress in M. japonicus, and that the hydrolysis of carbohydrates is increased as an adaptation to osmotic regulation during the early stages of stress. In addition, we found that phosphoenolpyruvate carboxykinase expression was increased significantly at 6, 12, and 24 h, after which it returned to normal, indicating that phosphoenolpyruvate carboxykinase is involved in the early stage of low-salinity stress in M. japonicus.

4.3 Immunity under low-salinity stress

Low-salinity stress can reduce the immunity of M. japonicus and increase its sensitivity to viral infection (Yu et al., 2003). KEGG analysis showed that the DEGs identified at 48 h were enriched in the Chagas disease (American trypanosomiasis) (ko05142) and Vibrio cholerae infection (ko05110) pathways, and most of these DEGs were downregulated, suggesting that M. japonicus is more resistant to pathogens under low-salinity stress. Shrimp lack an acquired immune system, and its innate immune system is divided into humoral immunity and cellular immunity. Shrimp are thought to be totally dependent on a congenital and nonadaptive mechanism to resist pathogen invasion (Li and Xiang, 2013; Tassanakajon et al., 2013; Wang et al., 2014a). KEGG analysis showed that the DEGs were enriched in lysosome (ko04142), peroxisome (ko04146), phagosome (ko04145), and apoptosis (ko04210) pathways related to cellular immunity, consistent with the changes in pathways identified previously in the L. vannamei hepatopancreas under low-salinity stress (Xu et al., 2017). Immune-related genes such as calnexin involved in the phagosome (ko04145) and cathepsin A involved in the lysosome (ko04142) showed significantly down-regulated expression in response to low-salinity stress, suggesting that low-salinity stress reduced phagocytosis and oxygen-independent reactivity in phagocytes of L. vannamei (Zhao et al., 2015). In this study, cathepsin C and lysosomal alpha-mannosidase in the lysosome (ko04142) and calreticulin in the phagosome (ko04145) were significantly downregulated. Of these genes, cathepsin C returned to normal after 12 h of down-regulation, lysosomal alpha-mannosidase returned to normal after 6 hours of down-regulation, and calreticulin was down-regulated only during the late stage (after 48 h) of stress, indicating that these genes are involved in different stages of cellular immunity and reduce the phagocytic capacity of shrimp. Apoptosis is a genetically programmed cell suicide process that eliminates unwanted or diseased cells and plays important roles in embryogenesis, homeostasis, insect metamorphosis and immunity (Steller and Grether, 1994; Opferman and Korsmeyer, 2003). Cytochrome C is a proapoptotic factor that is released from the outer surface of the inner mitochondrial membrane during the early stage of apoptosis. In combination with certain cytosolic proteins, it induces conversion of the latent apoptosis-promoting protease pro-caspase 9 into its active form (Skulachev, 1998; Marsden et al., 2002). In this study, we found that cytochrome C, involved in the apoptosis (ko04210) signaling pathway, was upregulated at 6 h; this is consistent with its previously reported up-regulation in L. vannamei infected with WSSV, indicating that cytochrome C may be induced by low-salinity stress to participate in apoptosis (Hu and Yao, 2016). Interestingly, we found that caspase 3 was down-regulated at 6 h, suggesting that M. japonicus responds to low-salinity stress by inhibiting caspase 3 and thus preventing apoptosis.

The humoral response includes the prophenoloxidase (proPO) system, hemolymph coagulation, and a variety of antimicrobial peptides (Li and Xiang, 2013; Tassanakajon et al., 2013; Wang et al., 2014a). Melanization is regulated by phenoloxidase (PO), which is activated by the proPO cascade control and plays an important role in the immune system of invertebrates, responding quickly to pathogen infection (Amparyup et al., 2013). Prophenoloxidase b has been shown to contribute to PO activity in crustacean hemolymph plasma (Masuda et al., 2012). We found that prophenoloxidase b was increased significantly at 12 and 24 h in oxidoreductase activity (GO: 0016491), suggesting that prophenoloxidase b promotes humoral immune regulation in M. japonicus under low-salinity stress. proPO activation and Toll pathway initiation in arthropods are mediated by serine proteinase cascades and regulated by serpins in the hemolymph (Liu et al., 2009). Injecting rPmSERPIN3 along with Vibrio harveyi into shrimp decreased the clearance rate of bacteria from the hemolymph, suggesting that PmSERPIN3 functions as a regulator of the proPO activation system (Wetsaphan et al., 2013). We found that the expression of serpin 3 was up-regulated after 6 hours of stress. We speculate that serpin 3 promotes activation of proPO and the Toll pathway during the early stage of low-salinity stress. In addition, antilipopolysaccharide factor 2, an antimicrobial peptide, was down-regulated in 6 and 12 h, indicating that anti-lipopolysaccharide factor 2 is also involved in immune regulation.

4.4 Chitin metabolic process under salt stress

Notably, the level of GO enrichment of DEGs related to chitin metabolism (q value: 6.4593E-10) was the highest among all metabolic processes. This indicates that chitin metabolism may be involved in osmotic regulation or affected by salinity stress, supporting the results of a salinity stress-based transcription analysis conducted in Portunus trituberculatus (Lv et al., 2013). In addition, we found DEGS in the control and low-salinity groups about the chitinase 1, peritrophin and peritrophin-like protein genes. In total, seven, five, and three chitinases were identified in L. vannamei, Penaeus monodon and M. japonicus, respectively. This study showed agreement that chitinase 1 is mainly expressed in the hepatopancreas (Watanabe et al., 1998; Huang et al., 2010; Proespraiwong et al., 2010; Rocha et al., 2012), and consistent with the data from this experiment. The growth, molting and development of crustaceans require chitinase, which functions in digesting chitinous food, modifying the parenteral nutrition membrane and degrading the chitin exoskeleton (Zou and Bonvillain, 2004; Zhang et al., 2014). In this study, chitinase 1 was initially down-regulated. We speculate that chitinase 1 is involved in the response to low-salinity stress and improves the ability of M. japonicus to resist low-salinity stress. The specific regulatory mechanism of this process remains to be explored. We presume that expression of chitinase elicits an immune response quickly under low-salinity stress, which promotes ecdysis, and then allows gradual adaption to changes in salinity, which may be related to osmotic regulation. Its subsequent inhibition under low-salinity conditions may be a regulatory mechanism to prevent energy waste. Peritrophin is the most widely investigated soluble protein in the insect peritrophic membrane. It has strong interactions with chitin-binding domain 2 and the envelope protein A domain to form the peritrophic membrane, which promotes food digestion and protects insects from microbial invasion (Jasrapuria et al., 2010). In recent years, many peritrophin-like genes have been identified and cloned in crustaceans such as F. merguiensis (Loongyai et al., 2007), P. monodon (Chen et al., 2009), F. chinensis (Du et al., 2006), and Exopalaemon carinicauda (Wang et al., 2013). Du et al. (2006) and Wang et al. (2013) found that silencing peritrophinlike genes by dsRNA interference improved the survival rate of shrimp infected with WSSV, suggesting that peritrophin-like genes may have an immune defense function. In our transcriptome data, chitinase 1 was initially up-regulated, then returned to normal, and then was down-regulated, while the peritrophin gene was initially down-regulated and then returned to normal, under low salinity compared with the control group. Interestingly, peritrophin-like protein gene showed the opposite trend of the peritrophin gene. We speculate that these three genes are involved in the regulation of low-salinity adaptation. Salt stress promotes a rapid immune response in shrimp. Shrimp may adapt to the low-salinity environment by selfregulating osmotic pressure. Down-regulation of chitinase levels may be undertaken to prevent energy waste. The specific regulatory mechanism of this response remains to be explored.

5 CONCLUSION

The transcriptome libraries in hepatopancreas of M. japonicus under low salt stress were successfully constructed, and many differentially expressed genes identified. The DEGs were annotated to obtain the pathways and genes related to osmoregulation, metabolism, and immunity. These functional genes provide valuable information for the further understanding of molecular mechanism of M. japonicus under low salt stress.

6 DATA AVAILABILITY STATEMENT

All data supporting the results of this study are available upon reasonable request from the corresponding author.