Transcriptome and Expression Profiling Analysis of the Hemocytes Reveals a Large Number of Immune-Related Genes in Mud Crab Scylla paramamosain during Vibrio parahaemolyticus Infection

Total RNA was extracted from the hemocytes of crabs at several time points, including 0, 12, 24, and 48 h after V. parahaemolyticus infection. Equal quantities of RNA were mixed together to construct a cDNA library for Illumina sequencing, which yielded a total of 62,004,374 reads. The reads containing more than five consecutive bases with a quality <13 were removed. The overall Illumina clean reads and bases are 52,934,042 and 4.76G, respectively. Files containing these data were deposited in the Sequence Read Archive of the National Center for Biotechnology Information (NCBI) with the accession number of SRP040563. The high-quality clean reads were assembled into 186,193 contigs by using the SOAP de novo software, which were further assembled into 81,709 consensus sequences with a mean length of 656 bp (range: 205∼23,885 bp). The length statistics of all consensus sequences is presented in Fig. 1.

59,120 unigenes were identified from 81,709 consensus sequences of mud crabs. After searching the Nr and Swissprot databases, 48,934 unigenes (31,482 in Nr and 27,498 in Swissprot) were identified. Among these, 10,566 unigenes belong to 3 categories of Gene Ontology, 25,349 to 30 categories of KEGG, and 15,191 to 25 categories of COG database, covering almost all functional categories (Table 1). The remaining 10,186 unigenes failed to match proteins in the Nr database or Swissprot databases and therefore represented potentially novel genes.

Gene ontology (GO) analysis was performed using the web-based Database for Annotation, Visualization, and Integrated Discovery (DAVID) [32]–[33]. Among the 59,120 assembled unigenes, 10,566 of them were successfully annotated by GO assignments to one or more of the three categories: biological process, cellular component and molecular function (Fig. 2), which fell further into 27, 17 and 17 subcategories with the largest ones being the “primary metabolic process” (59.1%), “cell” (53.9%) and “binding” (55.0%), respectively. In summary, these terms account for a large fraction of the overall assignments in the S. paramamosain transcriptome data.

Phylogenetic classifications of the predicted CDSs of unigenes were carried out by searching against the COG database (Fig. 3). A quantity of 15,191 unigenes were classified and subdivided into 25 COG categories (23 of them had known functions, while the other two were “general function prediction only” and “function unknown”), among which “General function prediction only” represented the largest group (9650 matched unigenes, 13.52%), followed by “Translation, ribosomal structure and biogenesis” (2,279, 14.08%) and “Cell wall/membrane/envelope biogenesis” (1,483, 9.16%).

A total of 25,349 unigenes in mud crab transcriptome were mapped to 30 statistically remarkable categories in KEGG (P<0.05) (S1 Table). Metabolic pathways (2701, 10.66%) and regulation of actin cytoskeleton (1610, 6.35%) were two of the most statistically significant categories. Notably, three immune related categories, mitogen-activated protein kinase (MAPK) signalling pathway (668, 2.64%), ubiquitin mediated proteolysis (594, 2.34%), and chemokine signalling pathway (761, 3.00%), were also identified as statistically significant.

To characterize the immune response of mud crab to bacterial infection, two DGE libraries were constructed using mRNA from the haemocytes of crabs infected with V.parahaemolyticus and the controls, respectively. After removal of low-quality tags, adaptor tags, and one copy number tag, a total of 7,562,126 and 7,468,808 clean tags were obtained from the libraries of bacterium infected and control crabs, respectively (S2 Table), which then were mapped to the transcriptome database described above. A total of 1213 differentially expressed genes (P<0.05) were found, including 538 up-regulated and 675 down-regulated genes, in the haemocytes of crabs infected with V.parahaemolyticus (Fig. 4). Files containing these differentially expressed genes upon V.parahaemolyticus infection were deposited in the Gene Expression Omnibus (GEO) of NCBI with the accession number of GSE45400↗.

To achieve a functional annotation of the infection-responsive genes, GO classifications were assigned to the differentially expressed genes by using DAVID, indicating that the up- and down-regulated genes were involved in immunity, transcription, translation regulations, and biological regulation (S3 Table). Among the 1213 differentially expressed genes, 961 (79.23%) genes were well annotated, whereas the remaining 252 (20.77%) genes had low sequence homology to known sequences in public databases, suggesting that they might be putative novel genes in S.paramamosain involved in immune response to V.parahaemolyticus infection.

Of these 1213 differentially expressed genes, 96 of them were known to be involved in immunity and immune related signal transduction in mud crabs after V. parahaemolyticus infection (S4 Table). Most of these immune-related genes were enriched in GO terms “Protein metabolic” and “immune system related”, which contained up to 40 and 20 genes, respectively. The rest three categories were “cytokine and chemokine” (18 genes), “Transporter and regulator activity” (10 genes), and “Cell proliferation and apoptosis” (8 genes).

Genes involved in the Toll pathway and the Immune deficiency (IMD) pathway, such as Toll-like receptor (TLR), interferon (IFN) and tumor necrosis factor receptor 6 (TNFR6), were found significantly induced in mud crab upon V. parahaemolyticus infection. Antimicrobial proteins, such as Crustin, ALF2, and ALF6, were also significantly up-regulated, indicating they may play important roles in crab's defence against V. parahaemolyticus infection.

Genes annotated as Ras superfamily, including Ras, Rab and myosin, were also activated by V. parahaemolyticus infection. Many genes in the transcription regulation group and genes encoding cell proliferation and apoptosis, including Tob1 and Syntaxin, as well as genes involved in antioxidant activity, such as AST, CAT, CathepsinA and CathepsinD, were found up-regulated in crabs infected with V. parahaemolyticus.

Of these differentially expressed genes, a large number of immune related genes belong to the proPO-activating system, in which the most significantly differentially expressed genes were serine proteinases, clip domain serine proteinase, serine protease-like protein, chymotrypsin, and their inhibitors, such as serpin and pacifastin.

Relative quantitative real-time RT-PCR (qRT-PCR) analysis was performed on randomly chosen genes (Chy, Fam, ALF2, Ras, Cd42, PpC, SS, Pero) to confirm the expression profiles of mud crabs in response to V. parahaemolyticus infection (Fig. 5). The results showed that qRT-PCR analyse of these genes was in good agreement with the Illuminia RNA-seq analyse, indicating the transcriptome data can reflect an actual gene expression profile in V. parahaemolyticus infected mud crabs.

Mud crabs used in this study were taken under permission from a local crab farm (Niutianyang, Shantou, Guangdong, China) (23.3543^oN. latitude and 116.6458^oE. longitude). Crabs were processed according to “the Regulations for the Administration of Affairs Concerning Experimental Animals” established by the Guangdong Provincial Department of Science and Technology on the Use and Care of Animals.

Samples collection and the pathogen challenge had been described previously [54]. Briefly, forty healthy mud crabs (100±10 g in weight) were collected from Niutianyang area in Shantou, China. Crabs were acclimated for 3 days in 1 m³ tanks (ten crabs per tank) in conditions (salinity: 8‰; temperature: 28°C) similar to those of culture ponds from which crabs were obtained. Sea water was changed twice a day. Crabs were then fed with shellfishes for one week before sampling. V. parahaemolyticus isolated from diseased crabs (from Shantou crab farms) were cultured in 2216E medium at 28°C. 30 crabs were inoculated individually in the ventral hind legs with 100 µL of V. parahaemolyticus (2.4×10⁷ cells/mL), while the other 10 crabs were injected with 100 µL of saline (0.9%) as the control group. 10 crabs were harvested for haemocyte sampling at 12 h, 24 h, 48 h after the challenge. Haemocytes were collected from blood samples after mixing with anti-coagulant solution and centrifugation at 800 g at 4°C for 20 min. Haemocytes were immediately frozen in liquid nitrogen and then stored at −80°C prior to RNA extraction.

Total RNA was extracted from 50 to 100 mg of tissue with TRIZOL Reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. The RNA samples were incubated for 30 min at 37°C with 10 units of DNase I (Takara, Dalian, China) to remove residual genomic DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260 nm/280 nm (A260/A280) using a ND-1000 spectrophotometer (LabTech, Holliston, MA) and integrity was ensured through analysis on a 1.5% (w/v) agarose gel.

cDNA libraries were constructed using RNA Sample Preparation Kit (Illumia, San Diego, CA) following manufacturer's instructions. Briefly, mRNA was purified from total RNA using polyATtract mRNA isolation systems (Promega, Madison, WI). 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 (Promega, Madison, WI). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H (Promega, Madison, WI). Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. In order to select cDNA fragments of preferentially 200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly). Products were quantified using the Agilent high sensitivity DNA assay on the Agilent Bioanalyzer 2100 system. Finally, the library preparations were sequenced on an Illumina Hiseq 2000 platform which generates 100 bp paired-end reads.

Raw sequencing reads were quality trimmed, and adaptor sequences were removed before the assembly. After removal of low quality reads, processed reads were assembled using SOAP de novo software [55] with default parameters. The overall assembly was performed using the combined sequence data for both treated crabs and the control groups. The contigs andsingletons were generally referred to as unigenes. As a result, 81,709 unigenes were generated. Subsequently, the unigenes were subjected to BLASTX similarity search against NCBI non-redundant protein database and the swissprot database using BLASTALL program with an E-value threshold of 10⁻⁵. All annotated unigenes were used to determine the COG term, GO term and KEGG pathway with a cut-off E-value of 10⁻⁵ using BLASTX [29]. The assembled contigs were used as a reference for annotating the DGE tags.

Gene expression was measured by counting tags from normal and pathogen-infected crabs and normalized to the total high-quality reads. High-throughput sequencing was performed using the Solexa/Illumina Genome Analyzer. The raw data was filtered to remove low quality sequences, including ambiguous nucleotides and adaptor sequences. The number of reads mapping to each unigene in the two populations analyzed (infected and non-infected) was counted and used as an approximate estimation of gene expression level in the corresponding individuals. For unigenes with multiple transcripts, the number of reads were summed and assumed to be the corresponding final read count. To link the expressed signatures to known genes from crab, the unigene dataset (http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?login=pda↗) from the transcriptome of S. paramamosain was used as a reference. The transcription level of each gene was deduced by determining the total number of reads mapped to each gene using Picard tools (http://picard.sourceforge.net/↗). To provide a relative assessment of transcript-abundance, the numbers of raw reads that were mapped to individual contigs were normalized for sequence length (i.e., Fragments Per Kilobase of transcript per Million mapped reads, FPKM). To investigate differences in gene expression profiles, FDR (false discovery rate) <0.001 was used as the threshold of P-value in multiple test to judge the significance of gene expression difference [56]. Genes were considered differentially expressed in a given library when the P-value <0.001 and a greater than two-fold change (absolute value of log2 ratio>1) in expression across libraries was observed.

Quantitative real-time PCR was performed using the ABI Prism 7300 Detection System (Applied Biosystems, Foster City, CA) with SYBR Green I as the fluorescent dye according to the manufacturer's protocol (Takara, Dalian, China). First-strand cDNA was synthesized from 2 µg of total RNA as described above and used as a template for real-time PCR with specific primers (S5 Table). Real-time PCR was performed in a total volume of 20 µl, and cycling conditions were 95°C for 5 min, followed by 40 cycles of 94°C for 5 s, 60°C for 31 s. All reactions were performed in biological triplicates, and the results were expressed relative to the expression levels of β-actin in each sample by using the 2^-ΔΔCT method [57]. Each sample was first normalized for the amount of template added by comparison with the abundance of β-actin mRNA [5], [58].

One way ANOVA analysis was used to analyze the means differences among samples. Pairwise t-test was then applied to compare gene expression levels between two crab populations. All statistical analyses (avo functions and pairwise t-test functions in R) were conducted in the R computation environment (http://www.r-project.org↗) [59]. P-values smaller than 0.05 were considered statistically significant.