Transcriptomic analysis of adaptive mechanisms in response to sudden salinity drop in the mud crab, Scylla paramamosain

Scylla paramamosain (Crustacea: Decapoda: Portunidae: Syclla De Hann) is a commercially important mud crab distributed along the coast of southern China and other Indo-Pacific countries (Lin Z, Hao M, Zhu D, et al, Comp Biochem Physiol B Biochem Mol Biol 208-209:29–37, 2017; Walton ME, Vay LL, Lebata JH, et al, Estuar Coast Shelf Sci 66(3–4):493–500, 2006; Wang Z, Sun B, Zhu F, Fish Shellfish Immunol 67:612–9, 2017). While S. paramamosain is a euryhaline species, a sudden drop in salinity induces a negative impact on growth, molting, and reproduction, and may even cause death. The mechanism of osmotic regulation of marine crustaceans has been recently under investigation. However, the mechanism of adapting to a sudden drop in salinity has not been reported. In this study, transcriptomics analysis was conducted on the gills of S. paramamosain to test its adaptive capabilities over 120 h with a sudden drop in salinity from 23 ‰ to 3 ‰. At the level of transcription, 135 DEGs (108 up-regulated and 27 down-regulated) annotated by NCBI non-redundant (nr) protein database were screened. GO analysis showed that the catalytic activity category showed the most participating genes in the 24 s-tier GO terms, indicating that intracellular metabolic activities in S. paramamosain were enhanced. Of the 164 mapped KEGG pathways, seven of the top 20 pathways were closely related to regulation of the Na+ / K+ -ATPase. Seven additional amino acid metabolism-related pathways were also found, along with other important signaling pathways. Ion transport and amino acid metabolism were key factors in regulating the salinity adaptation of S. paramamosain in addition to several important signaling pathways.


Background
Salinity as a key abiotic parameter that influences the distribution, abundance, physiology, and well-being of crustaceans. [3,4,15,34]. Salinity is also an important factor in the production of crustacean aquaculture [47], which can affect growth, survival, molting, oogenesis, embryogenesis and larval quality [6,14,21,30,32,35,38]. The salinity of crustacean aquaculture can easily be affected by a local torrential rain [46]. Fortunately, most marine species which have been studied have osmoregulatory capacities to adjust to the shifting salinities of estuary and wetland waters within limits. In crustaceans, the adaptability of osmoregulation is primarily achieved by the gills [7,12,16,28,34,36]. Since regulation of osmotic pressure in marine animals involves energy consumption, drastic changes in salinity can lead to death of the organism.
Recent studies have shown that low salinity influences ion channel activity [33,37,44,45] and L-type free amino acids [1,25,39,42,43], which are tightly involved with osmoregulation. In particular, the Na + / K + -ATPase, a well-known ion channel, is the main ion transport enzyme of post-larvae in crustaceans. Its function in the organism is to enhance adaptability to salinity changes through osmoregulation [20,22,26,27]. Chung and Lin [8] cloned the full-length α-subunit of the Na + / K + -ATPase cDNA, indicating the osmoregulatory role of the channel via both mRNA and protein expression. Likewise, Lu et al. [25] completed the cDNA cloning of glutamate dehydrogenase and its expression, indicating that GDH played an important role in controlling osmoregulation through free amino acids in S. paramamosain.
Although a great deal of progress had been made in the study of osmotic adjustment in crustaceans, there have been relatively fewer studies on S. paramamosain. In particular, the molecular mechanism of adaptation to sudden salinity drop has not yet been reported. Since S. paramamosain is a euryhaline species [40], it is easy to overlook the impact of salinity on the organism's physiology. In the environment, a sudden drop in salinity caused by heavy rain over a short period of time (drop by > 10‰) may lead to death. In this study, we simulated a drastic reduction in salinity from 23‰ to 3‰. Then, the molecular mechanism of adapting to the salinity drop was analyzed by transcriptome analysis.

Experimental animals and sectionalization
A total of 300 randomly selected crabs with a body weight of~30 g was selected and kept in a natural water environment with a salinity of 23 ‰ and a temperature of approximately 20°C. Every 50 crabs were randomly selected (weight~30 g) as a group, with a total of six groups, housed in six cement pools under identical physical and chemical conditions. The salinity of the seawater for three of the groups was adjusted to 3‰ from 23‰, which dropped by 20‰. These three groups were defined as the LS (low salinity) group. The other three groups were defined as the CK groups, where the salinity of seawater was kept at 23‰. All other conditions were the same as the LS group.
It should be noted that the space of a pond in the experiment was big enough for juvenile crabs, and we added a few tiles in ponds as shelter, which could effectively avoid fighting and killing each other.

HE staining
Gill morphology and ultrastructure of S. paramamosain were observed using light microscopy after hematoxylin (HE) staining. HE staining was conducted according to the method of Wang et al. [41]. First, gills were set in paraffin and sliced into sections with a thickness of up to 5-10 μm. Then, the sections were de-waxed using xylene and rehydrated in an ethanol series. Sections were stained with eosin and HE which purchased from Invitrogen (Carlsbad, CA, USA). Preparation of 4% paraformaldehyde solution were made using filtered isotonic sea water: the salinity of 23‰ and 3‰ filtered sea water was used in the C and LS groups, respectively.

Total RNA isolation and gene expression analysis
Total RNA was isolated from gill tissue using sqRT-PCR RNAiso Plus (TaKaRa, Dalian, China). The cDNA was synthesized using the Perfect Real Time version of the PrimerScriptTM RT reagent kit with gDNA Eraser (Perfect Real Time) (TaKaRa,) according to manufacturer's instructions. Then, sq-RT-PCR were chosen to analyze genes, and performed in a total reaction volume of 25 μl according to the manufacturers' instructions. The S. paramamosain beta-actin gene and 18S ribosomal RNA gene were selected as the internal control. Primers used in this study are listed in Additional file 1: Table S1.

Transcriptome sequencing
For Illumina paired-end sequencing, equivalent quantities of total RNA isolated from the three mud crabs were pooled as one sample, and eventually there were three samples in each group, CK and LS. After poly (A) mRNA was purified and fragmented into small pieces, we used random hexamer primers and reverse transcriptase (Invitrogen) to carry out first-strand cDNA synthesis. Second-strand cDNA synthesis was performed with RNase H (Invitrogen) and DNA polymerase I (New England BioLabs, Beijing, China). A cDNA library was constructed with average insert sizes of 200-500 bp and cDNA sequencing was conducted using the Illumina HiSeqTM 4000 system according to the manufacturer's protocols, with read length of 150 bp Transcriptome Quantification analyses two independent cDNA libraries were constructed for the two organs in parallel according to the Transcriptome protocol. The transcriptome sequencing was performed by BGI (BGI, Shenzhen, China).

Analysis of differentially expressed genes
Because the annotated genome of S. paramamosain has not been published, de novo assembly was used here as reference for further analysis. Firstly, the raw reads were filtered to remove adaptor and low quality sequences. Afterfiltration, clean reads were assembled into unigenes using Trinity de novo assembler, followed by TGICL clustering tool. The reads from control (CK) and experimental group (LS) were mapped against the assembled Unigene using HISAT. The FPKM (fragments per kb per million reads) method was used to calculate the expression abundance. Each unigene was subjected to a BLASTX search against the NCBI non-redundant (nr) protein database with an e-value threshold of 10-3. R package DESeq2 were performed to identify the differentially expressed genes (DEGs). DEG was considered as unigene with greater than 2-fold change and p-value < 0.05. Gene ontology (GO) terms and KEGG pathway annotation were achieved using the Blast2GO program and kaas (KEGG Automatic Annotation Server) on-line program (http://www.genome.jp/kaas-bin/kaas_main), respectively.

Results
Adaptive phenotype of S. paramamosain to sudden drop in salinity There were four deaths in the CK group within 7 days and 24 deaths in the LS group. The LS death time was concentrated within 24, 48, and 72 h. In addition, the LS group showed hyperactivity within 48 h. As time went by, the motility was diminished and normalized (Fig. 1a). The LS group did not have food over 72 h, and gradually started to eat over time. Conditions returned to normal after 120 h.
The gills are an important organ in osmoregulation of marine crustaceans, with the gill filament serving as the basic units of function. According to the results of gill slicing, under normal conditions, the gill filaments of the crab in the CK group were regular (Fig. 1b). Gills of crabs in the LS group became shorter and thicker from 6 h and reached the shortest at 72 h, which was less than half that of gills observed in the CK group ( Fig. 1a & b). Gills then gradually became longer with time and returned to normal after 120 h ( Fig. 1b & c). The changes in gill filament anatomy was consistent with those of the physiological activities mentioned above, suggesting that changes in gill filaments play an important regulatory role in adaptation of S. paramamosain to decrease in salinity.
Differentially expressed genes (DEGs) in the gill of S. paramamosain Gills of marine crustaceans play an important role in osmoregulation. In our study, the mud crabs had adjusted to a salinity of 3‰ in 120 h after a sudden drop in salinity. In order to study the molecular mechanism underlying this adaptation, we performed transcriptional profiling at 120 h. Approximately 39.6 Gb bases were generated in total on the BGISEQ-500 sequencing platform. Because the genome sequencing of S. paramamosain has not yet been elucidated, after reads filtering, Trinity [18] was used to perform de novo assembly with clean reads (Additional file 1: Table S2). Tgicl ( [31] Pertea et al., 2003) was used on cluster transcripts to remove abundance and get Unigenes (Additional file 1: Table S3). The proportion of bases with low quality (< 20) was very low in all samples, indicating a high quality of sequences. Finally, high-quality transcripts were obtained (Table 1) and used as reference sequences. Genes were annotated using Unigenes by aligning with seven functional database as follows 36 Figure S2). For functional annotation results, we detected 32,627 CDS by Transdecoder. We also detected 74,041 SSR distributed on 44,075 unigenes, and predicted 12,623 transcription factor (TF) coding unigenes.
A total of 249 genes was differentially expressed in the LS Group and the CK group, including 207 up-regulated genes and 42 down-regulated genes (fold change> = 2.00 Fig. 1 Food intake and morphology change of gills in response to the sudden drop in salinity in S. paramamosain. a, Food intake change in Scylla paramamosain. The juvenile crab with a body weight of~30 g eat about two Sinonovacula constricta (Lamarck) in 24 h, so three Sinonovacula constricta (Lamarck) per crab were feed every day and cleaned up the leftovers regularly to prevent the water going bad. b, the phenotypic changes of gill filaments within 120 h after the salinity drop: 23‰ in CK group, others were the phenotypes of 12, 24, 48, 72, 96, and 120 h after salinity dropped to 3‰. c, changes of gills in length; "0" indicated CK group, and the others were the gill length at 12, 24, 48, 72, 96, 120 h after salinity dropped to 3‰. Bars are 300 μm in b and adjusted p value <= 0.05) (Fig. 2). Of the 249 DEGs, 217 were annotated at least one of the following: Nr (192), Nt (160), Swissprot (154), KEGG (163), KOG (152), Interpro (130) and GO (32. No description was found for 3 DEGs (3 of which are down-regulated genes). Of the 192 differential genes annotated in the Nr database, 57 were described as "hypothetical" in the Gen-Bank database and were also excluded from further analysis. Finally, 135 DEGs that were annotated with Nr were screened, out of which 108 were up-regulated and 27 were down-regulated genes (Additional file 1: Table S4).

Functional annotation
Gene Ontology (GO) analysis was performed on DEGs using Blast2GO [9,17]. According to the second-tier GO terms, 32 DEGs were assigned 24 GO annotations which represented three main GO categories: biological process (18), cellular component (18), and molecular function (23) (Fig. 4 & Additional file 1: Figure S3). Among them, 24 GO terms were assigned to the up-regulated group and 17 to the down-regulated group. Twelve processes were identified in the biological process category, with 14 DEGs involved in cellular processes, 12 involved in metabolic processes, and 11 involved in single-organism processes. This identified three biological processes as the most strongly affected in the gill of S. paramamosain by the sudden drop in salinity (Fig. 4). In the cellular component category, membrane (14), membrane part (10), cell parts (7), and cell (7) were most involved (Fig. 4). In the molecular  (3) and electron carrier activity (1) (Fig. 4). It is worth noting that the catalytic activity category (14) had the most participating genes in the 24 s-tier GO terms. The results suggested that the salinity of the water environment dropped sharply. In order to maintain life activity, the metabolic activity of the cells in S. paramamosain was strengthened, in particular the catalytic function of enzymes regulating ion changes and osmotic pressure, processes involved in adapting to the new environment with dropped salinity.
In addition, among all the 24 s-tier GO terms, localization, response to stimulus, and multi-organization processes from the biological process category, as well as the extracellular region, membrane-enclosed lumen from cellular component, transporter activity, and electron carrier activity from molecular function categories were all upregulated DEGs with no downregulated DEGs.

Validity of DEGs in transcriptomic data
Ten differentially expressed genes were randomly sampled for verification by transcriptional level experiments, and the results were consistent with those of DEG analyses (Table 4 & Additional file 1: Figure S3), indicating that the results of DEG were reliable.

Discussion
The mechanism of osmotic regulation of aquatic crustaceans has received attention in recent years. Extensive work and quite a few important results have been reported on the morphological structure of osmotic regulatory organs [5,24], ion transport regulation [10,23,33], regulation of hemolymph osmoregulation ( [19];Huong et al.,   [42]), and neuroendocrine regulation [13,29]. For both euryhaline and stenohaline species, changes in salinity will result in the organism adapting through the regulation of the neuroendocrine system, osmotic regulatory organs (mainly gills), hematopoietic osmotic pressure, and ion transport. A series of changes will occur to adapt to the changing external environment in order to maintain normal physiological and metabolic activity. However, few studies on neuroendocrine regulation, regulation of ion transport enzymes, or osmotic regulation are reported for S. paramamosain. Thus, this paper focuses on osmotic regulation in aquatic crustaceans. Since sequencing of the genome of S. paramamosain has not been completed yet, much of the information obtained still depends on transcriptomics.
S. paramamosain is a euryhaline species, and especially loves living in shallow sea and estuary nearshore. In China, the salinity of seawater in S. paramamosain ponds on the farm is between 25‰ and 3‰ in most area, and the minority such as in Shanghai is below 3‰. There is a production experience in the actual production process that amplitude of variation in salinity exceeding 10‰ would cause death for mud crabs (This only refers to sudden salinity drop). 23‰ is a normal salinity of the seawater for juvenile crabs S. paramamosain, which was living in the salinity before our treatment. We had made a preliminary experiment to select a salinity for treatment with 10 juvenile crabs as a group, including 13‰, 8‰, 5‰, 3‰, and 1‰, and the degree of salinity drop was 10‰, 15‰, 18‰, 20‰, and 22‰, respectively. In the end, we found some individuals in salinity 3‰ begin to die, and most individuals nearly died in salinity 1‰. So 3‰ might be the optimal choice as a critical point for the research of adaptive mechanism responding to sudden drop in salinity, and we finally selected the 3‰ in our study. Moreover, the LS death time was concentrated within 24, 48, and 72 h. The LS group showed hyperactivity within 48 h, and as time went by, the motility was diminished and normalized. The LS group did not have food over 72 h, and gradually started to eat over time. Conditions returned to normal after 120 h. To be honest, the mechanism of adaptive process is very complex, and need more in-depth research in the future work. The study aimed at the adaptive mechanism in response to sudden salinity drop. So we compared the CK with 120 h_group (a state of complete adaptation) to discover the difference of between normal condition and the adaptive status after sudden salinity drop from 23‰ to 3‰, and eventually reveal the adaptive mechanisms in response to sudden salinity drop in the mud crab, S. paramamosain at the level of transcription.
To date, reports on the regulation of infiltration of S. paramamosain are still quite limited to cloning and expression of the Na + / K + -ATPase [8] and cloning and expression of glutamate dehydrogenase (GDH) [25]. GHD is an important enzyme for the metabolism of glycine, proline, and alanine, which serve as general osmolytes in aquatic animals. In this study, GO annotation of DEG analysis showed that the most involved genes were derived from the category of catalytic activity of molecular function (Fig.  4), suggesting that the possible role of free amino acids in osmotic regulation is by enzymolysis. In addition, we found that Na + / K + -ATPase is very active through the KEGG pathway classification and functional enrichment of DEGs (Figs. 5 and 6). The results showed that the Na + / K + -ATPase strengthened the ion exchange function necessary to maintain the osmotic balance required for normal survival after salinity drop from 23‰ to 3‰. In addition, there were many other KEGG pathways and differentially expressed genes, which might directly or indirectly participate in the regulation of osmotic adjustment of S. paramamosain, providing a valuable data source for subsequent studies.
S. paramamosain normally lives in estuary areas, where the environment is significantly different from freshwater, brackish water and seawater. Sometimes the salinity changes constantly, and this change requires a response in behavior, morphology, and biochemical physiology. In production, the sudden drop in salinity usually results from heavy rainfall in strong convective weather. For example, two typhoons hit Zhuhai in August 2017 in one week in Guangdong Province, causing huge losses to the crab farming industry. In this study, the cumulative salvage rate in six days was 16% (24/150) in the salinity sag test but may be higher in actual production. Because of the more complicated water environment system in production, the sudden drop in salinity is often accompanied by a decrease in water temperature. The comprehensive factors led to a decrease in crab immunity. As a result, there was an increase in pathogenic microorganisms in the water and an increase in crab mortality. In this study, only the salinity was changed and the rest of the environmental factors were controlled. For the first time, the molecular mechanism of S. paramamosain adapting to the salinity drop was studied. Through the research presented here, we had discovered a large number of potential genes that are related to the salinity adaptation in S. paramamosain. The possible KEGG pathway provided a basis for further research. In addition, this study was an important supplement to the physiological study of aquatic crustacean infiltration, but also provided a scientific basis for the regulation of crustacean aquaculture.

Conclusions
In conclusion, we analyzed transcriptomic changes in the gills after a sudden drop in salinity in S. paramamosain. One hundred thirty-five DEGs annotated by Nr were screened, of which 108 were up-regulated and 27 were down-regulated. GO analysis showed that catalytic activity (14) had the most participating genes in the 24 s-tier GO terms, indicating that intracellular metabolic activities in S. paramamosain were enhanced. Based on KEGG pathway and biological functional enrichment on DEGs, the top 20 pathways showed that seven pathways were directly related to the active regulation of the Na + /K + ATP enzyme: Proximal tubule bicarbonate reclamation (ko04964), protein digestion and absorption (ko04974), bile secretion (ko04976), thyroid hormone signaling pathway (ko04919), mineral absorption (ko04978), Insulin secretion (ko04911), thyroid hormone synthesis (ko0491), and the regulatory genes in Na + /K + ATPase were all up-regulated. Additionally, several amino acid metabolism pathways were detected: arginine biosynthesis (ko00220), alanine, aspartate and glutamate metabolism (ko00250), lysine degradation (ko00310), valine, leucine and isoleucine degradation (ko00280), amino sugar and nucleotide sugar metabolism (ko00520), biosynthesis of amino acids (ko01230), pyrimidine metabolism (ko00240). In addition, some famous signal pathways were found, such as cAMP signaling pathway (ko04024), MAPK signaling pathway (ko04013, ko04010), Wnt signaling pathway (ko04310), mTOR signaling pathway (ko04150), Ras signaling pathway (ko04014). Our findings suggest that not only Na + /K + ATPase and amino acids played a key role in osmoregulation, but also some important signal pathways participated in osmoregulation. Ultimately, survival of S. paramamosain may be sustained in new surroundings with a sudden drop in salinity (23‰ fell to 3‰), by adjusting to the low salinity. The functional genomic studies of DEGs obtained in this study allow for a better understanding of various physiological responses in marine crustaceans induced by a sudden drop in salinity.

Additional file
Additional file 1: Table S1. The gene-specific primers used in this study, Table S2. Clean reads quality metrics from the gill of S. paramamosain, Table S3. Quality metrics of unigenes from the gill of S. paramamosain ,  Table S4. DEGs annotation, Table S5. Human Diseases pathways and DEGs involved, Table S6. Environmental Information Processing pathways and DEGs involved, Table S7. Genetic Information Processing pathways and DEGs involved, Table S8. Metabolism pathways and DEGs involved, Figure S1. Distribution of base quality on clean reads from the gill of S. paramamosain, Figure S2. Venn diagram between NR, KOG, KEGG, Swissprot and Interpro. Figure S3. The distribution of DEGs in GO analysis, Figure S4. Validity of DEGs in Transcriptomic data. (DOCX 1505 kb) Abbreviations GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; LS: Low salinity; S. paramamosain: Scylla paramamosain; sqRT-PCR: semiquantitative reverse-transcription PCR