Comparative transcriptome analysis of Eriocheir japonica sinensis response to environmental salinity

Chinese mitten crabs (Eriocheir japonica sinensis) are catadromous, spending most of their lives in fresh water, but moving to a mixed salt-fresh water environment for reproduction. The characteristics of this life history might imply a rapidly evolutionary transition model for adaptation to marine from freshwater habitats. In this study, transcriptome-wide identification and differential expression on Chinese mitten crab groups were analysed. Results showed: clean reads that were obtained totalled 93,833,096 (47,440,998 in Group EF, the reference, and 46,392,098 in Group ES, the experimental) and 14.08G (7.12G in Group EF 6.96G in Group ES); there were 11,667 unigenes (15.29%) annotated, and they were located to 230 known KEGG pathways in five major categories; in differential expression analysis, most of the top 20 up-regulated pathways were connected to the immune system, disease, and signal transduction, while most of the top 20 down-regulated pathways were related to the metabolism system; meanwhile, 8 representative osmoregulation-related genes (14-3-3 epsilon, Cu2+ transport ATPase, Na+/K+ ATPase, Ca2+ transporting ATPase, V-ATPase subunit A, Putative arsenite-translocating ATPase, and Cation transport ATPase, Na+/K+ symporter) showed up-regulation, and 1 osmoregulation-related gene (V-ATPase subunit H) showed down-regulation. V-ATPase subunit H was very sensitive to the transition of habitats. These results were consistent with the tests of qRT-PCR. The present study has provided a foundation to further understand the molecular mechanism in response to salinity changing in water.


Introduction
Crustaceans are one of the most abundant and diverse animal groups, and they occupy diverse habitats [1][2][3][4]. It is well-known that some species such as the portunid live their whole lives in high-salt marine habitats. Interestingly, in the late Cretaceous, a marine crustacean ancestor a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 left the marine environment, transitioned to freshwater or terrestrial habitats, and then quickly occupied most of the world's rivers, lakes, rivulets, etc [5]. Now, some of them have completely adapted to their new low-salt environment and live in freshwater habitats throughout their lives (e.g., Freshwater crabs), while some of them are more specialized, spending most of their lives in freshwater, but moving to a mixed salt-fresh water environment for reproduction and larval development [6]. For the latter, in the reproductive season, they migrate to brackish estuary water to complete their habitat-mixing-dependent reproduction. The characteristics of this life history might imply a rapidly evolutionary transition model from freshwater to marine.
Adaptation to low-salt concentrations encountered by marine crustaceans was one of the most important challenges in their evolutionary history of transition from a marine to a freshwater environment. They must control osmoregulation to ensure the dynamic balance of Na + -K + inside and outside tissue cells [7][8]. Some shellfish species (e.g., Mytilus edulis, Scrobicularia plana, Mercenaria mercenaria) maintain the salinity balance by changing both the closing frequency of their double shell and the expression of transport enzymes [9]. To prevent unrestricted water diffusion, freshwater fishes expel excess water by producing a large amount of low-concentration urine through high-rate glomerular filtration, and reabsorb glucose and inorganic salt through proximal and distal tubules; the salt lost by urination will be regained through the gills and food [10]. In recent years, gene expression and regulation in several osmoregulation-related biomarkers (Na + /K + -ATP, 14-3-3, HSPs, NKCCla, Apo-14) were studied with the development of molecular techniques [11][12][13][14]. Among these, the Na + -K + -ATP enzyme was an important type of heterogenous transmembrane protease system involved in regulating the osmotic balance of aquatic crustaceans by improving enzyme activity and expression [13]. Generally, its molecular structure, expression, morphology and even salinity change regulation have been extensively studied [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. However, osmoregulation in aquatic crustaceans is a complex process with multi-proteins and multi-molecules. Furthermore, how to regulate the salt balance to rapidly adapt to increasing ambient salinity water habitats from fresh water was also a mystery for crustaceans (e.g., Chinese mitten crabs). Although, it had been reported that plasma cortisol concentration and gill Na + , K + -ATPase activity increased after exposure to high ambient salinity may be hypo-osmoregulatory mechanisms for whitefish (Coregonus lavaretus) [30], the osmoregulatory mechanism is very much species-specific [31]. Osmoregulation by Chinese mitten crabs (Eriocheir japonica sinensis) during physiological adaptation has been studied extensively, and recently genes (CA (carbonic anhydrase), CYP4C (cytochrome P450 4C), GDH (glutamate dehydrogenase), and the NHE (Na + /H + exchanger)) involved in osmoregulation had been also reported [31]. However, the genetic basis of osmoregulatory mechanisms in Chinese mitten crabs is still poorly understood.
Nowadays, analysis at the transcriptome level has become increasingly more popular to reveal molecular mechanisms of physiological change in organisms [32][33][34][35][36][37][38][39][40]. In this study, we first aimed to reveal several functional genes related to osmoregulation by transcriptome analysis through Illumina paired-end sequencing, assembly, annotation and expression of differential genes during the transition of freshwater and saline habitats. Then, we would provide more information about candidate genes associated with the osmoregulation response to salinity change in habitats of crustaceans.

Ethics statement
There are no specific permits for crab collection in the selected locations. The sampling locations are not privately owned or natural protected areas. Crabs used for the experiments are not considered endangered or protected species, and their collection is legal in China.

Animals collection and challenge experiments
In this study, Twenty-four adult individuals of E. j. sinensis with carapace length more than 55 mm [41] were collected from the aquatic market in September 2016 and acclimatized at 20˚C for 72 hours. Then, they were divided into two groups for construction of the comparative transcriptome database: EF (12) and ES (12). Group EF was cultured in freshwater for 120 hours, while group ES was challenged by a continuous increase of their water salinity: 5‰ (24 hours), 10‰ (24 hours), 15‰ (24 hours), 20‰ (24 hours) and 25‰ (24 hours). Finally, the three posterior gill tissues of the two groups were collected and frozen in liquid nitrogen for later experiments [42][43].

Library construction and Illumina sequencing
The total RNA was extracted from the gill mixture of every group using TRIzol 1 reagent (Invitrogen Corp., Carlsbad, CA, USA) according to the manufacturer's instructions. Quality of degradation and contamination was monitored on 1% agarose gels, and its purity was checked using the NanoPhotometer 1 spectrophotometer (IMPLEN, CA, USA). Concentration was measured using Qubit 1 RNA Assay Kit in Qubit 1 2.0 Fluorometer (Life Technologies, CA, USA); 1.5 μg RNA per sample was used as input material for the RNA preparations from gill samples. The integrity of RNA was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were generated using NEBNext 1 Ultra™ RNA Library Prep Kit for Illumina 1 (NEB, USA) following the manufacturer's recommendations and sequenced on an Illumina Hiseq 2500 platform with paired-end reads.

De novo transcriptome assembly
Raw reads were processed through Perl scripts, and then, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. At the same time, parameters such as Q20 and the GC-content of the clean data were calculated. Transcriptome assembly was accomplished using the Trinity programme (http://trinityrnaseq. sourceforge.net/) [44]. We obtained the transcripts by hierarchical clustering through Corset [45]. Finally, the longest transcripts for each gene were selected as unigenes. The value of N50 was also presented in this study to preliminarily estimate assembly quality.

Differentially expression analysis
The DESeq R package was used for differential expression analysis of two conditions/groups [47]. The p values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR) [48]. Genes with an adjusted p-value <0.05 found by DESeq were assigned as differentially expressed. GOseq R packages were implemented in GO enrichment analysis of the differentially expressed genes (DEGs)-based Wallenius non-central hypergeometric distribution. In addition, the KOBAS [49] software was used to test the enrichment significance of DEGs in KEGG pathways [50].

Expression analysis by quantitative real-time reverse transcription PCR (qRT-PCR)
Representative DEGs were selected for confirmation of RNA-Seq (Quantification) data by qRT-PCR using a SYBR Premix Ex Taq kit (Aidlab) according to the manufacturer's instructions. The 18S rRNA gene was used as an internal reference. Primers for qPCR were designed using Primer Premier 5.0 software. Reaction mixtures (20 uL) contained 10 μL of 2×SYBR 1 Premix Ex TaqTM buffer, 1 μL of forward and reverse primers, 1 μL of cDNA, and 7 μL of RNase-free H 2 O. The thermal profile for SYBR Green RT-PCR was 95˚C for 2 minutes, followed by 40 cycles of 95˚C for 15 seconds, 55-60˚C for 15 seconds, and 72˚C for 30 seconds. The relative gene expression was analysed using the comparative threshold cycle method [51].

Transcriptome sequence and assembly
Two cDNA libraries were constructed using gills from the freshwater group (EF) and the seawater group (ES) to obtain the transcriptome expression profile. Approximately 99,296,168 (49,845,580 in EF and 49,450,588 in ES) paired-end raw reads with a Q20% higher than 94.82% with a GC% of 42.15% and 46.02% in EF and ES, respectively, were obtained; This outcome could indicate good sequencing quality [52]. In this research, ultimately, 93,833,096 (47,440,998 in EF and 46,392,098 in ES) and 14.08G (7.12G in EF and 6.96G in ES) clean reads were obtained (File A in S1 File). The de novo transcriptome assembly by Trinity showed the total length and number of transcripts were 98,905,725 base pairs (bp) and 147,596 bp, respectively. The maximal length of transcripts was 18,715 bp (average length 670 bp; N50:1169 bp). The total length and number of unigenes were 79,287,774 bp and 76,307 bp, respectively (File B in S1 File). The maximal length of unigene was 18,715 bp with an average length of 670 bp. The parameter N50 for unigenes in this research was 1668 bp, which was slightly higher than that in other research also performed on crustaceans [53][54][55]. In this research, assembled sequences ranging from 300 to 1000 bp in length were the most abundant of the total sequence. A total of 10027 unigenes exceeded 2000 bp. For transcripts, the most abundant were clustered into a group of fewer than 1000 bp in length (File C in S1 File). In total, the assembly results also indicated that the length distribution pattern and mean length of the unigenes were consistent with previous transcriptome studies using Illumina Hiseq platform [56][57][58].

Unigene functional annotation and classification
Functional annotation. BLASTX was used to search the protein database, and 76,307 assembled unigene sequences were identified after transcriptome assembly (147,596  3,879 of these unigenes were annotated in all databases, which account for 5.08% of the total unigenes (File D in S1 File).
Transport and metabolism of ions should be important in osmoregulation, hence, unigenes that belong to the cluster of transport and metabolism were identified, including lipid transport and metabolism (505 unigenes), carbohydrate transport and metabolism (494 unigenes), and inorganic ion transport and metabolism (304 unigenes). Unigenes such as Na + / K + -ATPase, alpha subunit (Cluster-21379.0), multifunctional chaperone (14-3-3 family) (Cluster-173.0), aquaporin (major intrinsic protein family)(Cluster-18185.24311), Ca 2+ /Na + exchanger NCX1 and related proteins (Cluster-18185.34230) might be involved in the molecular osmoregulation of aquatic crustacean. These results provide an important and valuable database for future research into the molecular osmoregulation mechanism of aquatic crustacean acclimated to different habitats on salinity.

Differential expression and enrichment analysis
In KEGG analysis, the top 20 pathways were significantly involved in response to salinity changes (p<0.05) ( Table 1). Six pathways (30%) and a total of 102 unigenes were associated with metabolic processes (e.g., citrate cycle (TCA cycle), arachidonic acid metabolism, pancreatic secretion, glutathione metabolism, serotonergic synapse and nitrogen metabolism). Eight pathways (40%) and a total of 211 unigenes were associated with immune or disease processes (e.g., chemokine signalling pathway, B cell receptor signalling pathway, Fc gamma R-mediated phagocytosis, insulin signalling pathway, insulin resistance, platelet activation, shigellosis and bacterial invasion of epithelial cells). Thus, we have reason to believe that the unigenes in these pathways may play an important role in adapting to new environmental conditions during a habitat transition. Additionally, the top 20 up-regulated and top 20 down-regulated pathways of metabolic and immune systems were isolated. As a result, we found that most of the top 20 up-regulated pathways were related to the immune system, disease, and signal transduction, but nearly all the top 20 down-regulated pathways were related to metabolism systems ( Table 2). We inferred that pathways related to the immune system, disease, signal transduction, and metabolism should be given more attention for revealing mechanisms of rapid adaptation during the transition from low-salinity freshwater to marine habitats. Studies in the Comparative transcriptome analysis of Chinese mitten crab response to salinity short-term had partially revealed a relationship between metabolism and immunity that showed salinity (overnutrition) had impacts on crustaceans' immune defence abilities with a significant increase in the activity of immune enzymes [29]. Metabolism and immunity, undoubtedly, are inextricably linked and coordinated by common regulatory axes [59].

Osmoregulation-related genes
An interesting problem is discovering how osmoregulation-related genes are involved in ion transport and salt balance regulation in the transition from freshwater to saline water. According to the above results, most of the top 20 down-regulated pathways of Chinese mitten crabs were classified into metabolism pathways (some may relate to osmoregulation). However, most of the representative osmoregulation-related genes that encode the relative proteins (14-3-3 epsilon, Cu 2+ transport ATPase, Na + /K + ATPase, Ca 2+ transporting ATPase, V-ATPase subunit A, Putative arsenite-translocating ATPase, Cation transport ATPase and Na + /K + symporter) were up-regulated (p<0.05). Only the V-ATPase subunit H was down-regulated (p<0.05) ( Table 3). The results of real-time RT-PCR showed that expression of DEGs annotated as Cluster-18185.13163 (14-3-3 epsilon), Cluster-18185.30584 (Cu 2+ transport ATPase), Cluster-18185.15624 (Na + /K + ATPase, beta subunit), Cluster-18185.21625 (Ca 2+ transporting ATPase), Cluster-18185.25046 (Na + /K + ATPase, alpha subunit), Cluster-18185.26540 (V-ATPase subunit A), Cluster-18185.27366 (V-ATPase subunit D), Cluster-18185.8525 (Ptype ATPase), Cluster-18185.921 (Clchannel proteins) were up-regulated in treatment group compared to the control. Osmoregulation-related genes in the process of habitat transition included both up-regulation and down-regulation. These results were consistent with those of Hui et al. (2014), who performed a desalination-based transcriptome analysis on larvae of Chinese mitten crabs. Membrane permeability is related to osmoregulation [60]; for this reason, up-regulation of 14-3-3 epsilon protein may have an important role in osmoregulation because the 14-3-3 epsilon protein could suppress activation of Bcl-2, an antagonist of cell death, and further stop initiating mitochondrial membrane permeability [61]. Moreover, transporters are important osmoregulation factors responsible for the uptake and efflux of important substances such as inorganic ions, sugars and amino acids. This study also reflected that the Cu 2+ transport ATPase, Na + /K + ATPase, beta subunit, Ca 2+ transporting ATPase, Cation transport ATPase, etc. were down-regulated in this study when adult Chinese mitten crabs were moved . This outcome implies sensitivity to ion flow regulation as part of V-ATPase to adaptation. Overall, this rapid adaptation to habitat transition must be related to many more pathways, such as environmental information processing, signal transduction, and energy metabolism.
In the late Cretaceous period, an ancestor of marine crustaceans left the long-term safety of the marine environment, transitioned to freshwater or terrestrial environments, and occupied most of the world's rivers, lakes, and rivulets, etc. During their evolutionary history, crustacean ancestors regulated the balance of Na + -K + inside and outside of the tissue cells by gene overexpression or mutation. In the present study, several functional genes related to osmoregulation have been revealed using transcriptome analysis and testing by qRT-PCR (Fig 4). We showed  these genes might have played an important role in habitat transition of crustaceans. Perhaps, molecular adaptive evolution in the long term happened in some osmoregulation-related genes in the form of an increasing number of subunits, and these subunits then played an important role in rapid adaptation to the freshwater for Chinese mitten crabs. As transcriptome information from more crustacean species is accumulated, more detail in their rapid adaptation and molecular mechanisms will be available in the future.
Supporting information S1 File. Statistics of transcriptomic profile on quality, assembly, length distribution and annotation. (RAR)