Transcriptome Analysis of Portunus trituberculatus in Response to Salinity Stress Provides Insights into the Molecular Basis of Osmoregulation

Background The swimming crab, Portunus trituberculatus, which is naturally distributed in the coastal waters of Asia-Pacific countries, is an important farmed species in China. Salinity is one of the most important abiotic factors that influence not only the distribution and abundance of crustaceans, it is also an important factor for artificial propagation of the crab. To better understand the interaction between salinity stress and osmoregulation, we performed a transcriptome analysis in the gills of Portunus trituberculatus challenged with salinity stress, using the Illumina Deep Sequencing technology. Results We obtained 27,696,835, 28,268,353 and 33,901,271 qualified Illumina read pairs from low salinity challenged (LC), non-challenged (NC), and high salinity challenged (HC) Portunus trituberculatus cDNA libraries, respectively. The overall de novo assembly of cDNA sequence data generated 94,511 unigenes, with an average length of 644 bp. Comparative genomic analysis revealed that 1,705 genes differentially expressed in salinity stress compared to the controls, including 615 and 1,516 unigenes in NC vs LC and NC vs HC respectively. GO functional enrichment analysis results showed some differentially expressed genes were involved in crucial processes related to osmoregulation, such as ion transport processes, amino acid metabolism and synthesis processes, proteolysis process and chitin metabolic process. Conclusion This work represents the first report of the utilization of the next generation sequencing techniques for transcriptome analysis in Portunus trituberculatus and provides valuable information on salinity adaptation mechanism. Results reveal a substantial number of genes modified by salinity stress and a few important salinity acclimation pathways, which will serve as an invaluable resource for revealing the molecular basis of osmoregulation in Portunus trituberculatus. In addition, the most comprehensive sequences of transcripts reported in this study provide a rich source for identification of novel genes in the crab.


Introduction
Portunus trituberculatus (Crustacea: Decapoda: Brachyura), commonly known as the swimming crab, is widely distributed in the coastal waters of Korea, Japan, China, and southeast Asia [1]. This species inhabits estuaries and coastal waters, which belong to typical euryhaline crab species. In China, it is a major edible crab species and one of the most important fishery resources [2] and the production has now reached 90,000 tons, valued at AUS$2.5 billion in 2009 [3]. salinity is one of the most important abiotic factors that influence not only the distribution and abundance of crustaceans, but also their general physiology and well being [4].The water salinity condition is also an important factor for artificial propagation of the swimming crab [5]. Throughout their prolonged culture period, Portunus trituberculatus often experience substantial salinity fluctuations either due to heavy rainfalls or droughts, which could have significant impacts to farm productivity and in severe situations, heavy mortality. This often requires the ability to regulate hemolymph osmolytes in relation to the environment they inhabit via osmoregulation to control their hemolymph osmotic pressure [6].
Due to the implications and critical importance of osmoregulation to the crab artificial propagation, a number of researchers have been devoted to this topic. An extensive literature that describes the growth, development, physiology, behavior, and propagation techniques of Portunus trituberculatus exposed to salinity stress have revealed the crab grows in optimal salinity ranged from 20-35ppt, whereas they can occur at salinities below 6 ppt and will survive salinities in excess of 48ppt [7][8][9][10][11][12]. In order to study the mechanism of osmoregulation, Xu et.al. investigated gene expression in the Portunus trituberculatus exposed to different salinity stresses via cDNA microarray chip, and 417 differentially expressed genes were identified [5]. Their study revealed a few important salinity acclimation pathways, which may be helpful in understanding the molecular basis of osmoregulation and salinity adaptation in the crab.
Even though cDNA microarray technology is a powerful tool for studying genome-wide gene expression, this technology fails to detect sequence variation and to recognize new genes or transcripts and can only be designed from limited expressed sequence tag data as the genome of Portunus trituberculatus has not yet been determined. To date, there are 13,985 ESTs available for the crab in the Genebank, however, it remains insufficient for the comprehensive understanding of Portunus trituberculatus transcriptome. Many low expression transcripts would be missed from current EST data, which makes it difficult for further analysis on transcriptome.
Next-generation high-throughput RNA sequencing technology (RNA-seq) is a recently-developed method for discovering, profiling, and quantifying RNA transcripts with several advantages over other expression profiling technologies including higher sensitivity and the ability to detect splicing isoforms and somatic mutations [13]. Because it is not restricted by the unavailability of a genome reference sequence, this approach has been applied in decoding the genomes of several non-model organisms, providing valuable information in the understanding of gene function, cell responses and evolution [14][15][16]. Significant progress has also been made in understanding the transcript expression of various marine crustacea by RNA-seq over the last two years, such as Litopenaeus vannamei, Fenneropenaeus chinensis, Eriocheir sinensis and Macrobrachium nipponense [17][18][19][20][21].The countable, almost digital, nature of RNA-seq data makes them particularly attractive for the quantitative analysis of transcript expression levels, which can give reliable measurements of transcript levels in one or more conditions [22]. However, such investigations in Portunus trituberculatus have not been reported.
In the present study, we examined the whole transcriptome responses to salinity stress of the Portunus trituberculatus for the first time using the Illumina's sequencing technology. Considering individual monitoring of the Portunus trituberculatus responses to salinity stress, nine libraries (three technical replicates of each condition) were established from the gills of Portunus trituberculatus that exposed to optimal, low and high salinity seawater, respectively. The study aimed to compare the expression patterns of the three conditions to better understand the transcriptomic regulation in Portunus trituberculatus to salinity stress and identify genes involved in osmoregulation of the crab. The results of this study are an important resource for future researches on mechanism of osmoregulation for marine invertebrates.

Salinity Challenge Experiment and Sample Preparation
The swimming crabs, Portunus trituberculatus at 80 days age (5.62~11.66g in body weight),were obtained from a local farm in Qingdao, China. All the samples were acclimated in the laboratory (33 ppt, 18°C) for one week before the experiment treatment. The crabs were divided into 3 groups (90 crabs for one group) and acclimated to low salinity challenged (LC, 5 ppt), non-challenged (NC, 33 ppt), high salinity challenged (HC, 50 ppt) at 18°C. For each group, the sixth pair of gills from 9 crabs were collected after ten days and samples were stored in RNAlater (Ambion) at 4 °C over night and then at -20 °C until RNA extraction within 2 weeks.

RNA Isolation, cDNA Library Construction and Illumina Deep Sequencing
Total RNA was isolated from each sample by trizol (Invitrogen, CA,USA). RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the Nano Photometer® spectrophotometer (IMPLEN, CA, USA) . RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA).RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).A total amount of 3g RNA per sample was used as input material for the RNA sample preparations and all samples had RIN values above 8. Then, samples of three individuals were pooled within each group in equal amounts to generate three mixed sample. The pooling samples were then used to prepare 9 separate Illumina sequencing libraries (three technical replicates of each condition). cDNA libraries were generated using Illumina TruSeq™ RNA Sample Preparation Kit (Illumia, San Diego, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in Illumina proprietary fragmentation buffer. First strand cDNA was synthesized using random oligonucleotides and SuperScript II. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities and enzymes were removed. After adenylation of 3' ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 200 bp in length the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 10 cycle PCR reaction. Products were purified (AMPure XP system) and quantified using the Agilent high sensitivity DNA assay on the Agilent Bioanalyzer 2100 system.
In the final step before sequencing, all 9 individual libraries were normalised and pooled together in a single lane on an Illumina HiSeq2000 platform and 90 bp paired-end reads were generated.

Bioinformatic Analysis
Quality control. Raw data (raw reads) of fastq format were firstly processed through our self-written perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GCcontent and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.
Transcriptome assembly. The left files (read1 files) from all libraries/samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file.Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity [23] with min_kmer_cov set to 2 and all other parameters set default.
Transcriptome Annotation. The unigenes were compared with the protein nonredundant database using BlastX with E values less than 1.0×10 -5 (E values less than 1.0×10 -5 were considered as significant) [24]. All annotated unigenes were used to determine the Clusters of Orthologous Groups of proteins (COG) term, Gene Ontology (GO) term and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway with a cut-off E-value of 1.0×10 -5 using BlastX.
Differential Expression, Cluster analysis and GO enrichment Analysis. Differential expression analysis was performed using the DEGseq (2010) R package. P value was adjusted using q value [25]. qvalue<0.005&|log2 (foldchange)| >1was set as the threshold for significantly differential expression. The three samples from each treatment were used to generate mean expression levels. Hierarchical cluster analysis is used to identify differentially expressed genes with certain patterns of expression under two different salinity challenge using Cluster 3.0 [5,26]. The settings for the calculations were as follows [5]: similarity is measured by standard correlation, and clustering method is average linkage. The genes are clustered according to similarities in expression profiles across three salinity conditions. GO enrichment analysis of the differentially expressed genes (DEGs) was performed using GOseq based Wallenius non-central hypergeometric distribution [27], which can adjust for gene length bias in DEGs.

Real-time RT-PCR confirmation of Illumina sequencing data
To validate our Illumina sequencing data, twelve differentially expressed genes were selected for quantitative RT-PCR analysis, using the same RNA samples as for the transcriptome profiling. Primers were designed using the Primer5 software (Premier Biosoft International) (Table S1).
RpL8 gene was selected as an internal control for qPCR analysis and the primers reference Xu's literature [5]. First strand cDNA was synthesized from1 mg of RNA using M-MuLV reverse transcriptase (Qiagen). The qPCR reaction mixture (20 mL) consisted of 26 Power SYBR Green PCR Master mix, 0.9M each of the forward and reverse primers, and 1 mL of template cDNA. PCR amplification was performed under the following conditions: 50°C for 2 min and 95°C for 30 s, followed by 40 cycles of 95°C for 15 s and 62°C for 1 min, and a final extension at 72°C for 5 min.

Result Illumina Draft Reads and Sequence Assembly
In this study, gills of three crabs were used to prepare one pooled RNA sample for each group of LC, NC and HC (three technical replicates of each group). Nine cDNA libraries were then constructed to perform Illumina sequencing. The schematic of Illumina deep sequencing and analysis are shown in Figure 1. We obtained 27,696,835, 28,268,353, and 33,901,271 qualified Illumina read pairs for LC, NC and HC, giving rise to total clean bases of 5.60, 5.72 and 6.80G, respectively. The overall Illumina read pairs and clean bases for all samples are 89,866,459 and 18.16G, respectively ( Table  1). Files containing these data were deposited in the Short Read Archive of the National Center for Biotechnology Information (NCBI) with accession numbers of SRR1013694 (LC), SRR1013695 (NC) and SRR1013696 (HC).
After assembly analysis based on all Illumina reads, we identified 141,339 contigs. The average length of all contigs was 1,157 bp, with the smallest sequence of 201 bp and the largest one of 36,385 bp. The sequence length distribution of contigs is indicated in Figure 2 and Table 2. Further assembly analysis showed that all contigs contributed to 94,511 unigenes, with an average length of 644 bp.
To assess the abundance and coverage of the transcriptome data, we matched the assembled unigenes against the known EST library on Genbank. The 13,985 ESTs downloaded from NCBI were clustered and assembled, and 2,612 assembled EST-unigenes with mean length of 783 bp were generated. Comparisons between transcriptome unigenes and ESTunigenes were performed using BLASTn algorithm. Results were shown in Figure 3 as a Venn chart. 74.3% (1,942 out of 2,612) of the EST-unigenes can be matched in the transcriptome unigenes library, whereas only 2.05% (1,942 out of 94,511) of the transcriptome unigene sequences can be found in the ESTs library.

Annotation of Unigenes
After ruling out short-length and low-quality sequences, 94,511 unigenes were selected and subjected to annotation analysis by matching sequences against Nr and Swiss-prot databases using BLASTx searching with an E value 1.0×10 -5 and all annotated unigenes were used to determine the COG term, GO term and KEGG pathway ( Table 3 and Figure S1, S2, S3). 13,212 unigenes (13.98% of the total) can be matched in Nr database, and 11,711 (12.39% of the total) matched in Swissprot. For main species distribution matched against Nr

Differentially Expressed Genes
Comparison of gene expression showed that a total of 1,705 unigenes were differentially expressed after two different salt challenges (qvalue<0.005 & |log2 (foldchange)|>1), including 615 (158 up and 457 down) differentially expressed unigenes between NC and LC (Table S2), and 1516 (895 up and 621 down) between NC and HC (Table S3). Moreover, 426 unigenes were significantly differentially expressed in both NC Vs. LC and NC Vs. HC (Figure 4). To validate our Illumina sequencing results, tweleve differentially regulated genes were selected for quantitative real time-PCR (qRT-PCR) analysis, of which ten genes agrees well with the illumina sequencing analysis.The other two, including sodium-potassium-chloride cotransporter and V-type proton ATPase subunit F, do not match the Illumian sequencing analysis perfectly, however, it shows similar trends in up-or down-regulation of genes analysised by Illumina sequencing (table 4). The results suggesting that the majority of the salinity responsive genes found from the Illumina sequencing data analysis are authentic.
The differentially expressed genes were further categorized into eight patterns based on expression profiles via hierarchical cluster analysis ( Figure 5): genes highly up-regulated in low and high salinity conditions (Cluster I, 39 genes), only in low salinity challenge (Cluster V, 104 genes) and only in high salinity conditions (Cluster VII, 840 genes). These categories are in contrast with three clusters that showed different kinetics of down-regulation in the two salinity conditions (Cluster II, 356 genes), only in low salinity challenge (Cluster VI, 85 genes), and only in high-salinity environment (Cluster VIII, 250 genes). Moreover, Cluster III (15 genes) showed the genes highly upregulated in low salinity and down-regulated in high salinity challenge, while Cluster IV (16 genes) exhibited the genes highly up-regulated in high salinity and down-regulated in low salinity challenge.

Discussion
In recent years, the next generation sequencing methods have also been applied to analyze transcriptomes of crustaceans using Illumina sequencing. Comparing with the traditional methods, the next-generation high-throughput DNA sequencing techniques provide more ideal methods for transcriptome analyses with high efficiency, low cost and high data output. The former studies on Portunus trituberculatus transcriptome were performed using traditional cDNA library and Sanger sequencing methods with RNA from many organs such as gill, eyestalk, blood and hepatopancreas [28][29][30][31], however, it remains insufficient for the comprehensive understanding of Portunus trituberculatus transcriptome. In our study, we generated 141,339 contigs of Portunus trituberculatus transcriptome based on the next generation sequencing techniques. Further assembly analysis showed that all contigs contributed to 94,511 unigenes, largely enriching the transcriptome data of Portunus trituberculatus and prompting the genome studies of crustaceans.
We compared our transcriptome data with Portunus trituberculatus EST sequences obtained from NCBI and showed that more than half of the EST sequences (74.3%) can be matched in the transcriptome data, whereas only 2.05% of the transcriptome unigene sequences can be found in the ESTs library. It suggests the transcriptome data provide abundant information besides the now available ESTs sequences. After a homology search in the non-redundant protein database at NCBI, a total of 13,212 unigenes, which took up a proportion of 13.98% in all the unigenes, showed significant BlastX hits of known protein sequences. The distribution of significant BlastX hits over different organisms was also analyzed. Due to the lack of genomic information in Portunus trituberculatus, the majority of the assembled sequences (11.25%) matched genes from microcrustacean arthropod Daphnia pulex which have full genomic information.
Comparison of gene expression among the different treatment groups in the current experiment is helpful for identification of candidate genes underlying response to salinity stress in Portunus trituberculatus. In this study, we detected a total of 1,705 unigenes were differentially expressed between the comparison of LC vs NC and HC vs NC (qvalue<0.005 & | log2 (foldchange)|>1), which significantly more abundant than the data of previous study, in which only 417 differentially expressed genes were found after two different salt challenges (10% and 40%) via cDNA microarray technology [5]. The results prove that the next generation sequencing methods should be more powerful than cDNA microarray technology in expression analysis because it can provide high data output and recognize new unigenes or unique isoforms present in transcriptome [32].
Many fewer differentially expressed unigenes were found in the LC Vs. NC (615) compared with that of HC Vs. NC (1,516). To further unravel the significantly altered biological processes upon salinity stress, the up-and down-regulated genes categorized into eight patterns based on expression profiles were subjected to the GO term enrichment analysis. The results revealed that a total of 454 unigenes were enrichmented in 63 processes, among which, the differentially expressed unigenes in low salinity stress were enrichmented in 25 processes (6 processes were induced, 17 processes were suppressed and 2 processes were were affected). The differentially expressed unigenes in high salinity stress were enrichmented in 58 processes (31 processes were induced, 24 processes were suppressed and 3 processes were were affected). Further analysis found, 5 and 38 processes were only enrichmented in low salinity stress or high salinity stress respectively, another 20 processes were enrichmented in both stress. In conclusion, our analysis uncovered: 1, more genes and processes were involved in high salinity stress adaptation than low salinity stress adaptation. 2, the processed involved in low salinity stress were mainly supressed, and the processed involved in high salinity stress were mainly induced. 3, hypoosmoregulation and hyper-osmoregulation mechanisms share some processes but also has many differences.
Past research has revealed two major strategies, i.e. the "limiting process" and the "compensatory process", that are adopted by crustaceans for osmoregulation and both are predominately accomplished by the gills [33,34]. A "limiting process"is a strategy whereby the maintenance of hemolymph osmolality/ions is accomplished by adjusting the permeability of the boundary structures (e.g. gill mem-branes), which can be a highly effective method to reduce ion diffusion and water inux rather than solely relying on the more energetically demanding mechanisms of ion transport. It has been suggested that the mechanisms during long term salinity exposure is often the result of gill membrane fatty acid compositional changes [35]. In the case, generally membranes containing higher proportions of fatty acids with high unsaturation indexes (i.e. higher double bonds). In this research, we found at least five up-regulation unigenes belong to fatty acid biosynthetic process (GO:0006631) in response to salinity stress (Table S4) which suggested the strategies of "limiting process" maybe play an important role in osmoregulation in Portunus trituberculatus.
A "compensatory process", on the other hand, is a strategy whereby the maintenance of hemolymph osmolality/ions is accomplished via the active movement of solutes into or out of the hemolymph to counter balance their passive diffusion [4,33]. When the environmental osmolality is higher than that of the hemolymph, it is termed "hyper-osmotic" and crustaceans must compensate for ion influx via hypo-osmoregulation; for the opposite, when the environmental osmolality is lower than that of the hemolymph, it is termed "hypo-osmotic" and crustaceans must compensate for ion loss from their hemolymph via hyper-osmoregulation [34]. Portunus trituberculatus belongs to weak hyperosmoregulators [36] which employs the apical Na +, K + , 2Cl( -) cotransporter driven by the inwardly directed Na + gradient, supplemented by apical K + channels, that recycle K + and hyperpolarize the apical membrane, creating a negative cell potential that drives Clefflux through basal Clchannels. Na + uptake may be Figure 5. Eight expression patterns ofthe differentially expressed genes. Cluster I, genes highly up-regulated in low and high salinity conditions; Cluster II, genes down-regulated in low and high salinity conditions;Cluster III, genes highly up-regulated in low salinity and down-regulated in high salinity challenge; Cluster IV, genes highly up-regulated in high salinity and down-regulated in low salinity challenge; Cluster V, genes highly up-regulatedonly in low salinity challenge; Cluster VI, genes down-regulation only in low salinity challenge; Cluster VII, genes highly up-regulatedonly in high salinity challenge; Cluster VIII, genes down-regulatedonly in high salinity challenge. augmented through the apical sodium/hydrogen exchanger. The resulting outside-positive transepithelial potential drives substantial inward paracellular Na + movement across this leaky epithelium [36,37]. In our study, ion transport-related genes including chloride channel protein, sodium/hydrogen exchanger, sodium/glucose cotransporter, and carbonic anhydrase were up-regulated in response to low salinity stress. In addition, we found that five signifantly up-regulated expressed unigenes were enrichmented in anion transport process (GO:0006820). The differential expression of this key ion transport genes potentially supports this models of decapod Crustacea hyper-osmoregulation reported previously [36,37].
On the contrary, in hypo-osmoregulation, passive NaCl influx is compensated by active NaCl secretion, however, the mechanisms of active NaCl secretion in the Crustacea are not known [36]. In this studys, we also found ion transport-related genes in high salinity stress, but most of them down-regulated expression which were different from low salinity stress. These down-regulated genes contained Na + , K + , 2Cl( -) cotransporter, carbonic anhydrase, calcium-activated chloride channel regulator 1, and anion exchange protein etc. Based on previous researches [36,37] and our results, we speculated that: in hypo-osmoregulation, for maintaining hemolymph osmotic balance, Portunus trituberculatus initiatively downreglated the expression of some important ion transport-related genes to prevent high concentration of ions of the external environment penetrate into the body. However, more detailed studies are required to disclose the mechanisms of hypoosmoregulation in Crustacea.
The aquaporins were considered as important salt responsive effectors [38]. The functions of AQP have been studied in mammals [39] and fish [40] which suggested that epithelial water flux occurs, at least in part, through transcellular pathways formed via aquaporins (AQPs), however, the functions of aquaporins in osmoregulation was little known in crustacea. Previous researches suggested that during acclimation to a hyperosmotic environment, increasing renal AQPs expression may play a role in enhancing water reabsorption by the kidney [41] and translational gene knockdown of AQPs protein reduced water influx in fish [42]. Additionally, the water channel activity of AQPs can be modulated through phosphorylation and dephosphorylation [43][44][45]. In addition to maintaining ionic homeostasis under salt stress conditions, Portunus trituberculatus also need establish water homeostasis via aquaporins (AQPs). Our data showed that the mRNA expression level of AQPs decreased under hypo-osmotic and hyper-osmotic stress conditions. Similar results were found in many previous studies under long-term saltstress [38,46]. These results may indicated that long-term stress induced a reduction of AQPs activities via transcriptional level controls and post-translational modifications to protect against water currents as well as gill cell swelling and shrinkage. However, AQPs activities whether affected by Figure 6. GO enrichment analysis of the differentially expressed genes. Cluster I, genes highly up-regulated in low and high salinity conditions; Cluster II, genes down-regulated in low and high salinity conditions;Cluster III, genes highly up-regulated in low salinity and down-regulated in high salinity challenge; Cluster IV, genes highly up-regulated in high salinity and down-regulated in low salinity challenge; Cluster V, genes highly up-regulatedonly in low salinity challenge; Cluster VI, genes down-regulation only in low salinity challenge; Cluster VII, genes highly up-regulatedonly in high salinity challenge; Cluster VIII, genes down-regulatedonly in high salinity challenge.  It has been reported that Na + ,K + -ATPase is the driving force in establishing an ion gradient across the epithelial cell membrane in marine crabs [47]. In this study, Na + ,K + -ATPase β-subunit showed significant down-regulation in high and low salinity conditions, which suggested that Na + ,K + -ATPase might play a very important role in salinity adaptation. Previous studies have shown that Na + ,K + -ATPase β-subunit showed significant up-regulation during low salinity challenge and significant down-regulation during high salinity challenges [5], which was different from our data. This may be due to different salinity challenged in the two studies. Although lots of literature showed that gene expression levels of Na + ,K + -ATPase α subunit were highly up-regulated during salinity stress [48][49][50], significant up-/down-regulation of Na + ,K + -ATPase α subunit was not detected in our transcriptome data, and similar results were also found in Xu's studies [5]. None of the significant gene expression changes of Na + ,K + -ATPase α subunit during salinity challenges do not negate its role for osmoregulation in swimming crab, and more future researches should be done toward Na + ,K + -ATPase α subunit.
Previous studies have shown that the gill apical V-H + -ATPase in freshwater crabs was shown to be involved in ion regulation, and marine ones generally showed a cytoplasmic V-H + -ATPase distribution, which did not participate in osmoregulation [51]. Similar suggestions have already been proposed in studies of the gills of Eriocheir sinensis and Uca tangeri [52][53][54]. Portunus trituberculatus belong to typical euryhaline marine crab species, however, it is interesting that four V-H+-ATPase subunits were showed significantly different expresstion in high salinity stress, which suggested V-H + -ATPase might play a very important role during hypoosmoregulation. Subcellular localization and function of V-H + -ATPase deserve further study in Portunus trituberculatus.
The present investigation suggests that the degradation of hemolymph or muscle protein to FAAs or de novo synthesis of free amino acids (FAA) may serve as a mechanism to compensate or alleviate the effects of ion influx under salinity stress [55], and that glycine, proline and alanine widely function as osmoeffectors in a number of crustacean species [56][57][58][59]. It must be noted, however, such a response is not universal to all crustaceans [60]. In this study, many differentially expressed unigenes were enrichmented in amino acid metabolism and synthesis processes. Among which, glycine catabolic process (GO:0006546) was induced in low salinity stress, and in high salinity stress, glycine metabolic process (GO:0006544) and Lserine metabolic process (GO:0006563) were suppressed, but lysine biosynthetic process (GO:0009089) was induced. Interestingly, we aslo found many differentially expressed unigenes enriched in proteolysis process (GO:0006508), which down-regulated (11 genes) in low salinity stress and mainly upregulated (43 up-regulated genes and 20 down-regulated genes) in high salinity stress. This results uncovered that, in high salinity stress, to balance the osmotic pressure inside and outside the body, Portunus trituberculatus increase levels of FAAs in vivo via the way of increasing protein degradation, decreaseing amino acid catabolism, and increasing synthesis of certain non-essential amino acids, and in low salinity stress, Portunus trituberculatus takes the opposite approach to hypertonic regulation. In addition to the known processes associated with osmoregulation in Crustacea., there were some other processes enriched which were little known in osmoregulation. The top three processes based on the number of unigenes were oxidation-reduction process (GO:0055114), cellular protein modification process (GO:0006464) and phosphorylation (GO:0016310), and the top three processes based on the enrichment level (p-value, excluding the processes of less than three unigenes) were chitin metabolic process (GO:0006030), carbohydrate metabolic process (GO: 0005975) and small GTPase mediated signal transduction (GO:0007264).
The oxidation-reduction process is essential in maintaining cellular homeostasis. Under physiologic conditions, cells maintain oxidation-reduction balance through generation and elimination of reactive oxygen/nitrogen species (ROS/RNS) [61]. Normally, the oxidation-reduction process homeostasis ensures that the cells respond properly to endogenous and exogenous stimuli. However, when the oxidation-reduction process homeostasis is disturbed, oxidative stress may lead to aberrant cell death and contribute to disease development [61]. Both exogenous and endogenous sources contribute to the formation of intracellular ROS/RNS. This reseach suggested that exposure to salinity stress, especially the high salinity stress has been shown to active the oxidation-reduction process in Portunus trituberculatus in order to balance the level of ROS/RNS. Of course, there was another possibility that the oxidation-reduction process in Portunus trituberculatus has been disturbed in this high salinity stress. Both cellular protein modification and phosphorylation process belong to the scope of post-translational modifications.
Posttranslational modifications modulate the activity of most eukaryote proteins [62] which. manifest as chemical modifications that occur on amino acid side chains in a sitespecific manner. They can temporarily or permanently change the fate of the protein by enhancing the functionality and/or stability of the target protein through the recruitment of auxiliary factors, change the proteins' cellular localisation or signal the most terminal fate, proteasomal degradation [63]. Our data showed that many gene were enrichmented in cellular protein modification (50 unigenes ) and phosphorylation process (31 unigenes ), suggested that posttranslational modifications play an important role in osmoregulation of Portunus trituberculatus, and further research is worthwhile.
Many differentially expressed unigenes (38) enriched in chitin metabolic process with the maximum enrichment level (pvalue: 1.60E-08) were aslo found in our data, which suggested chitin metabolic process may be involved in osmoregulation or affected by salinity stress. Further analysis found there were 4 chitinase genes. In crustaceans, during the molting cycle, chitinase dissolves chitin in the old exoskeleton into more soluble forms, which can then be partially reabsorbed into the body and used to synthesize the new exoskeleton [64]. Therefore, the chitinase are indispensable for crustaceans. By now, many genes encoding chitinases have been isolated and characterized in crustaceans [64][65][66][67] due to their importance in growth, development and immunity. However, chitinase with function of osmoregulation has not been reported in crustacea.
More detailed research can be done according to the connections between the salinity changes and chitin metabolic process in the future.
Small GTP-binding genes play a crucial regulatory role in a number of cellular processes in both plants and ani-mals, such as vesicle-mediated intracellular trafficking, signal transduction, cytoskeletal organization, and cell division [68]. Differential expression of small GTP-binding proteins has been previously reported in various abiotic stresses in plant [68], animal [69], and bacteria [70], however, the functions of Small GTP-binding genes was little known in crustacea. 15 unigenes were enriched in GTPase mediated signal transduction process with a relatively low p-value implied that GTPase mediated signal transduction process plays an important role in Portunus trituberculatus salinity adaptation. In addition, 17 up-regulated unigenes were enriched in carbohydrate metabolic process in high salinity stress suggested that hypo-osmoregulation of Portunus trituberculatus depends on energy consumption. This result have also been reported in Xu's studies [5] and further experiments are needed to verify the hypothesis.
It should be noted that, the processes of response to oxidative stress (GO:0006979), defense response (GO: 0006952), response to oxidative stress (GO:0006979) and defense response to bacterium (GO:0042742) were enriched in salinity stess, and lots of stress-related and immunity-related genes such as HSP, cathepsin, lectin, serine protease, and peroxiredoxin show significantly up-regulated or downregulation expression responding to low or high salinity stress, which suggested that salinity dilution or elevation appeared to invoke a classic ''stress response'' and "immunity response" at the transcriptional level in Portunus trituberculatus. This results similar to Xu et al.(2011)'s results [5] and different from Towle et al.(2011)'s results [71], which might reflect that Portunus trituberculatus possesses lower salinity adaptability than C. maenas, and it also suggested that the two crab species might have different signal pathways against the salinity ''stress'' [5].
It is noteworthy that the present study identified a large number of transcripts (526 unigenes) significantly upregulated or downregulate during salinity stress for which no annotation was readily available. This new genes provide a wealth of reference data to further research the mechanism of osmoregulation in Crustacea.

Conclusions
In general, our work represents the first report of the utilization of the next generation sequencing techniques for the study of osmoregulation in Portunus trituberculatus. Through more than 18.16G clean bases obtained, we established transcriptome data set for Portunus trituberculatus subjected to salinity stress, and the data we generated could enrich on genomic resources of this non-model organism. Based on the comparison of gene expression, a substantial number of genes were found to be modified by salinity stress which demonstrated the complexity of salinity adaptation mechanism in the crab. Our results revealed a few important salinity acclimation pathways via enrichment analysis, which may be helpful in understanding the molecular basis of osmoregulation in Portunus trituberculatus. Table S1. Primer sequences and product size of target and reference genes used for quantitative RT-PCR. (XLSX )   Table S2.

Supporting Information
List of differently expressed genes from Portunus trituberculatus in response to low salinity challenged (LC) versus non-challenged (NC).