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-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-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-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.

Materials and Methods

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, GC-content 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 hyper-geometric 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 ). 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 . 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 ). 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).

Figure 1

Schematic of Illumina deep sequencing and analysis.

It includes sample preparation, cDNA library construction and Illumina sequencing, data analysis including assemble, blast, GO annotation, gene expression analysis, etc.

Table 1

Summary of draft reads of three libraries by Illumina deep sequencing.

Groupa	PE library size (bp)	Read pairs	Read length (bp)	Cleanbases (G)
LC	200	27,696,835	90	5.60
NC	200	28,268,353	90	5.72
HC	200	33,901,271	90	6.84
Total	200	89,866,459	90	18.16

a, the data of every group comes from three technical replicates

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 and Table . Further assembly analysis showed that all contigs contributed to 94,511 unigenes, with an average length of 644 bp.

Figure 2

Sequence length distribution of contigs assembled from Illumina reads.

Table 2

Assembly Results Statistics.

Transcript length interval	200-500bp	500-1kbp	1k-2kbp	>2kbp	total
Number of transcripts	73,467	24,848	19,568	23,456	141,339

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 EST-unigenes were performed using BLASTn algorithm. Results were shown in Figure 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.

Figure 3

Venn chart for comparisons between assembled transcriptome unigenes and assembled EST-unigenes.

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 and Figure ). 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 database, 11.25% of the matched unigenes showed similarities with Daphnia pulex, followed by Tribolium castaneum (6.90%), Pediculus humanuscorporis (5.41%), Branchiostoma floridae (3.98%), Paramecium tetraurelia (3.58%), Strongylocentrotus purpuratus (3.42%), Nasonia vitripennis (3.17%), Ixodes scapularis (2.90%), Megachile rotundata (2.84%), Acyrthosiphon pisum (2.66%), Camponotus floridanus (2.11%), Harpegnathos saltator (2.07%), Danaus plexippus (1.98%), Pediculus humanus (1.57%), Tetrahymena thermophila (1.36%), Aedes aegypti (1.36%), and others (43.45%).

Table 3

Summary statistics of Portunus trituberculatus transcriptome blast assignment.

Unigenes	Nrannotations	Swiss-protannotations	ESTscan prediction	COGhits	GOmapped	KEGG hits
94,511	13,211	11,711	39,625	11,528	23,914	5,419

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 ), and 1516 (895 up and 621 down) between NC and HC (Table ). Moreover, 426 unigenes were significantly differentially expressed in both NC Vs. LC and NC Vs. HC (Figure ). 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 ). The results suggesting that the majority of the salinity responsive genes found from the Illumina sequencing data analysis are authentic.

Figure 4

Differentially expressed genes that are unique or shared among LC vs NC and HC vs NC.

Table 4

Gene expression levels of 12 selected genes for the comparison of controland salinity-challenged groups (5% or 45%)according to Illumina sequence and qPCR analysis.

Gene	ID	qPCR			Illumina sequence
		LC vs. NC	HC vs. NC		LC vs. NC	HC vs. NC
glutamine synthetase	comp33831_c0	0.81	0.20*		0.73	0.16*
glutamate dehydrogenase	comp40590_c0	1.10	0.31*		1.20	0.36*
Sodium / potassium ATPase beta chain	comp53011_c0	0.11*	0.27*		0.13*	0.25*
aquaporin	comp54248_c0	0.32*	0.10*		0.23*	0.01*
heat shock protein 70	comp54992_c0	0.53	4.26*		0.64	2.35*
carbonic anhydrase	comp55558_c0	7.36*	0.54*		4.73*	0.44*
V-type proton ATPase subunit B	comp62089_c0	1.11	0.39*		0.59	0.46*
sodium-potassium-chloride cotransporter	comp70314_c0	4.13*	0.31*		1.03	0.26*
chitinase	comp71131_c0	0.34*	0.16*		0.23*	0.04*
V-type proton ATPase subunit F	comp74516_c3	3.46*	5.39*		1.45	2.13*
Chloride channel protein	comp76355_c5	6.71*	1.63		2.80*	0.85
sodium/hydrogen exchanger	comp76376_c2	9.10*	1.35*		2.56*	1.81

Asterisk (*) marks the significant difference betweensalinity-challenged group and control group (P < 0.05).

The differentially expressed genes were further categorized into eight patterns based on expression profiles via hierarchical cluster analysis (Figure ): 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 up-regulated 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.

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.

GO enrichment analysis

Biological process of the differentially expressed genes of the eight patterns were further analyzed according to GO functional enrichment analysis (Figure , Table and Table ). GO functional enrichment analysis results showed cluster I were involved in chitin metabolic process (GO:0006030, 3 unigenes). Cluster II were involved in ten biological processes at a significant level. The mainly processes and involved unigene numbers are chitin metabolic process (GO:0006030, 13 unigenes), protein polymerization (GO:0051258, 6 unigenes), response to oxidative stress (GO:0006979, 4 unigenes), potassium ion transport (GO:0006813, 4 unigenes). Cluster III were involved in four biological processes which mainly include oxidation-reduction process (GO:0055114, 4 unigenes), lipid metabolic process (GO:0006629, 3 unigenes). Cluster IV were involved in five biological processes which mainly include oxidation-reduction process (GO:0055114, 5 unigenes), proteolysis (GO:0006508, 3 unigenes). Cluster V were involved in three biological processes which mainly include defense response (GO:0006952, 5 unigenes), anion transport (GO:0006820, 5 unigenes). Cluster VI were involved in five biological processes which mainly include proteolysis (GO:0006508, 8 unigenes), chitin metabolic process (GO:0006030, 3 unigenes). Additionally, cluster VII and cluster VIII were involved in twenty-nine and sixteen biological processes at a significant level respectively. The mainly processes, involved unigene numbers and cluster numbers are oxidation-reduction process (GO:0055114, 72 unigenes, cluster VII), cellular protein modification process (GO:0006464, 50 unigenes, cluster VII), proteolysis (GO:0006508, 40 unigenes, cluster VII), phosphorylation (GO:0016310, 31 unigenes, cluster VII), oxidation-reduction process (GO:0055114, 27 unigenes, cluster VIII), proteolysis (GO:0006508, 20 unigenes, cluster VIII), carbohydrate metabolic process (GO:0005975, 17 unigenes, cluster VIII), chitin metabolic process (GO:0006030, 10 unigenes, cluster VIII).

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.

Table 5

GO enrichment analysis of the differentially expressed genes.

Cluster	GOID	Term	Level of GO	Number of unigenes	P-value
I	GO:0006030	chitin metabolic process	6	3	0.001204934
II	GO:0006030	chitin metabolic process	6	13	1.60E-08
II	GO:0051258	protein polymerization	7	6	0.000898448
II	GO:0006813	potassium ion transport	8	4	0.032189626
II	GO:0006979	response to oxidative stress	3	4	0.000251565
II	GO:0006120	mitochondrial electron transport, NADH to ubiquinone	9	3	0.037277062
II	GO:0006875	cellular metal ion homeostasis	8	3	0.041235929
II	GO:0009755	hormone-mediated signaling pathway	7	3	0.037277062
II	GO:0016568	chromatin modification	6	3	0.043988317
II	GO:0071941	nitrogen cycle metabolic process	3	3	0.005597731
II	GO:0006935	chemotaxis	4	2	0.048209006
III	GO:0055114	oxidation-reduction process	3	4	0.017862676
III	GO:0006629	lipid metabolic process	3	3	0.035596858
III	GO:0006546	glycine catabolic process	10	2	9.13E-06
III	GO:0008033	tRNA processing	9	2	0.002337968
IV	GO:0055114	oxidation-reduction process	3	5	0.043166178
IV	GO:0006508	proteolysis	5	3	0.031323303
IV	GO:0002098	tRNA wobble uridine modification	12	2	0.000102853
IV	GO:0000154	rRNA modification	10	2	0.000190739
IV	GO:0009452	7-methylguanosine RNA capping	9	2	0.000190739
IV	GO:0001510	RNA methylation	8	2	0.000395991
IV	GO:0006479	protein methylation	8	2	0.000770984
V	GO:0006820	anion transport	6	5	5.52E-05
V	GO:0006952	defense response	3	5	0.034822636
V	GO:0006614	SRP-dependent cotranslational protein targeting to membrane	10	2	0.008743629
VI	GO:0006508	proteolysis	5	8	0.003397869
VI	GO:0006030	chitin metabolic process	6	3	0.01432079
VI	GO:0006888	ER to Golgi vesicle-mediated transport	7	2	0.037838935
VI	GO:0009878	nodule morphogenesis	6	2	0.021455881
VI	GO:0001539	ciliary or flagellar motility	5	2	0.001800326
VII	GO:0055114	oxidation-reduction process	3	72	0.000952872
VII	GO:0006464	cellular protein modification process	6	50	0.000230686
VII	GO:0006508	proteolysis	5	40	0.000397024
VII	GO:0016310	phosphorylation	5	31	0.001296757
VII	GO:0008610	lipid biosynthetic process	4	19	0.003426418
VII	GO:0009966	regulation of signal transduction	5	17	0.041422297
VII	GO:0007264	small GTPase mediated signal transduction	6	15	0.000113889
VII	GO:0006030	chitin metabolic process	6	11	0.006871104
VII	GO:0009108	coenzyme biosynthetic process	5	10	0.001612753
VII	GO:0006888	ER to Golgi vesicle-mediated transport	7	7	0.021852563
VII	GO:0006006	glucose metabolic process	7	6	0.038801546
VII	GO:0042398	cellular modified amino acid biosynthetic process	9	4	0.04770803
VII	GO:0035023	regulation of Rho protein signal transduction	8	4	0.008616499
VII	GO:0032446	protein modification by small protein conjugation	8	4	0.039101467
VII	GO:0007339	binding of sperm to zona pellucida	7	4	0.000990193
VII	GO:0006099	tricarboxylic acid cycle	7	3	0.047820393
VII	GO:0006000	fructose metabolic process	7	3	0.007238233
VII	GO:0001503	ossification	3	3	0.043900004
VII	GO:0009089	lysine biosynthetic process via diaminopimelate	11	2	0.04770803
VII	GO:0015813	L-glutamate transport	11	2	0.018737173
VII	GO:0006425	glutaminyl-tRNA aminoacylation	11	2	0.018737173
VII	GO:0009168	purine ribonucleoside monophosphate biosynthetic process	10	2	0.041097287
VII	GO:0034968	histone lysine methylation	10	2	0.03486768
VII	GO:0046058	cAMP metabolic process	10	2	0.02365927
VII	GO:0006570	tyrosine metabolic process	9	2	0.014309718
VII	GO:0006541	glutamine metabolic process	9	2	0.04770803
VII	GO:0006784	heme a biosynthetic process	8	2	0.010408479
VII	GO:0009231	riboflavin biosynthetic process	7	2	0.029045762
VII	GO:0030534	adult behavior	4	2	0.03486768
VIII	GO:0055114	oxidation-reduction process	3	27	0.001223294
VIII	GO:0006508	proteolysis	5	20	1.19E-05
VIII	GO:0005975	carbohydrate metabolic process	3	17	6.23E-05
VIII	GO:0006030	chitin metabolic process	6	10	2.79E-06
VIII	GO:0006544	glycine metabolic process	9	4	0.000286422
VIII	GO:0006631	fatty acid metabolic process	8	4	0.013582145
VIII	GO:0042742	defense response to bacterium	5	4	0.020716006
VIII	GO:0006825	copper ion transport	9	3	0.002850419
VIII	GO:0007218	neuropeptide signaling pathway	7	3	0.042788588
VIII	GO:0006979	response to oxidative stress	3	3	0.002600038
VIII	GO:0006334	nucleosome assembly	10	2	0.01916756
VIII	GO:0009396	folic acid-containing compound biosynthetic process	10	2	0.01916756
VIII	GO:0006563	L-serine metabolic process	9	2	0.00077586
VIII	GO:0009084	glutamine family amino acid biosynthetic process	9	2	0.005774778
VIII	GO:0008272	sulfate transport	8	2	0.011002118
VIII	GO:0007631	feeding behavior	3	2	0.020714011

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-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, hypo-osmoregulation 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 ) 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 Cl- efflux through basal Cl- channels. Na+ uptake may be 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 down-reglated 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 hypo-osmoregulation 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-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 phosphorylation in Portunus trituberculatus should be further verified.

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-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-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 hypo-osmoregulation. 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-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 L-serine 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 up-regulated (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 site-specific 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 (p-value: 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-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 down-regulation 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.

Primer sequences and product size of target and reference genes used for quantitative RT-PCR.

(XLSX)

Click here for additional data file.

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

(XLSX)

Click here for additional data file.

List of differently expressed genes from Portunus trituberculatus in response to high salinity challenged (HC) versus non-challenged (NC).

(XLSX)

Click here for additional data file.

COG Classification of the unigenes. Possible functions of 11528 unigenes were classified and subdivided into 26 COG categories.

(JPG)

Click here for additional data file.

GO classification of all unigenes. Most unigenes canbe divided into three major categories, including biological process, cellular component, and molecular function.

(JPG)

Click here for additional data file.

KEGG Classification of the unigenes. 5419 unigenes were assigned into 31 KEGG pathways.A, Cellular Processes; B, Enviromental Information Processing; C, Genetic Information Processing; D, Metabolism; E, Organismal Systems.

(PNG)

Click here for additional data file.