In-depth comparative transcriptome analysis of intestines of red swamp crayfish, Procambarus clarkii, infected with WSSV

Crayfish has become one of the most important farmed aquatic species in China due to its excellent disease resistance against bacteria and viruses. However, the antiviral mechanism of crayfish is still not very clear. In the present study, many high-quality sequence reads from crayfish intestine were obtained using Illumina-based transcriptome sequencing. For the normal group (GN), 44,600,142 high-quality clean reads were randomly assembled to produce 125,394 contigs. For the WSSV-challenged group (GW), 47,790,746 high-quality clean reads were randomly assembled to produce 148,983 contigs. After GO annotation, 39,482 unigenes were annotated into three ontologies: biological processes, cellular components, and molecular functions. In addition, 15,959 unigenes were mapped to 25 different COG categories. Moreover, 7,000 DEGs were screened out after a comparative analysis between the GN and GW samples, which were mapped into 250 KEGG pathways. Among these pathways, 36 were obviously changed (P-values < 0.05) and 28 pathways were extremely significantly changed (P-values < 0.01). Finally, five key DEGs involved in the JAK-STAT signaling pathway were chosen for qRT-PCR. The results showed that these five DEGs were obviously up-regulated at 36 h post WSSV infection in crayfish intestine. These results provide new insight into crayfish antiviral immunity mechanisms.

invading pathogens via an efficient and specific immune pattern 8 . The study of the intestinal transcriptomes is an important part of the research of the innate immune response mechanism.
In recent years, next-generation sequencing (NGS) technology has been widely used to screen out large amounts of genetic information in model organisms 9 . NGS technology is superior in many aspects to the traditional Sanger sequencing technology. NGS technology can provide enormous amounts of sequence data in a much shorter amount of time and at a much cheaper cost 10 . The expressed sequences produced using NGS technology are usually ten-fold or one-hundred-fold greater than those that are produced using traditional Sanger sequencing technology 11 . In the present study, Hi-Seq technology was used to sequence the intestinal transcriptomes of crayfish from a normal group (GN) and a WSSV-challenged group (GW). This information was used to generate expression profiles and to discover differentially expressed genes between normal crayfish and WSSV-challenged crayfish. Moreover, the functions of differently expressed genes (DEGs) were annotated and classified by the Gene Ontology (GO) database, Clusters of Orthologous Groups (COG) database, and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. These data provide an important resource for research on the gene functions, molecular events, and signaling pathways relating to the invertebrate antivirus immune response.

Results and Discussion
Illumina sequencing of the crayfish intestinal transcriptome. Illumina-based RNA sequencing was carried out with two sets of crayfish intestine samples (GN and GW). After cleaning and quality testing the GN sample, a total of 44,600,142 clean reads were screened out from 46,945,132 raw reads, corresponding to 4,460,014,200 total nucleotides (nt). The Q20 percentage (percentage of bases whose quality was greater than 20 in clean reads), N percentage (percentage of uncertain bases after filtering), and GC percentage were 98.05%, 0.00% and 40.45%, respectively (Table 1). For the GW sample, a total of 47,790,746 clean reads were screened out from 49,574,674 raw reads, corresponding to 4,779,074,600 nt. The Q20 percentage, N percentage, and GC percentage were 97.96%, 0.00% and 41.58%, respectively ( Table 1). All of these sequences were used for further analysis.
De novo assembly of the transcriptome. After low-quality reads and short reads were removed, high-quality clean reads were used to carry out transcriptome de novo assembly using Trinity software with the default parameters 12 . For the GN sample, a total of 44,600,142 high-quality clean reads were randomly assembled to produce 125,394 contigs with an N50 of 701 bp. The contigs were further assembled and clustered into 70,791 unigenes with a mean length of 627 nt that were composed of 7,352 distinct clusters and 63,439 distinct singletons. In addition, the N50 of the above unigenes was 1,258 bp ( Table 2). For the GW sample, a total of 47,790,746 high-quality clean reads were randomly assembled to produce 148,983 contigs with an N50 of 769 bp. The contigs were further assembled and clustered into 83,043 unigenes with a mean length of 672 nt that were composed of 9,686 distinct clusters and 73,357 distinct singletons. In addition, the N50 of the above unigenes was 1,456 bp ( Table 2). The distributions of the unigene sequence lengths for the GN and GW samples are shown in Figs 1 and 2, respectively.
Functional annotation of predicted proteins. After transcriptome de novo assembly for two sets of crayfish intestine samples (GN and GW), the transcripts were used for annotation, in combination with previously reported data from two other transcriptomes 13 . First, sequence annotation was carried out based on unigenes from the merged group 14 . Then, the putative functions of all of the unigenes were analyzed based on GO and COG classifications 15 . In this study, before the analysis of DEGs associated with white spot syndrome virus (WSSV) infection, a basic sequence analysis of the merged group transcriptome data (98,676 unigenes) was  In addition, a GO classification analysis was carried out; GO classification is an internationally standardized gene function classification system. This analysis provides a dynamically updated controlled vocabulary and can exactly define gene characteristics and products in any organism 16 . This analysis includes three ontologies: biological process, cellular component, and molecular function 17 . The biological process ontology includes biological adhesion, biological regulation, cell killing, cellular component organization or biogenesis, cellular  process, developmental process, establishment of localization, growth, immune system process, localization, locomotion, metabolic process, multi-organism process, multicellular organismal process, negative regulation of biological process, negative regulation of biological process, positive regulation of biological process, regulation of biological process, reproduction, reproductive process, response to stimulus, rhythmic process, signaling, and single-organism process. The cellular component ontology includes cell, cell junction, cell part, extracellular matrix, extracellular matrix part, extracellular region, extracellular region part, macromolecular complex, membrane, membrane part, membrane-enclosed lumen, organelle, organelle part, synapse, synapse part, virion, and virion part. The molecular function ontology includes antioxidant activity, binding, catalytic activity, electron carrier activity, enzyme regulator activity, molecular transducer activity, nucleic acid binding transcription factor activity, protein binding transcription factor activity, receptor activity, structural molecule activity, and transporter activity.
In the present study, a GO analysis was carried out using Blast2GO software. A total of 13,848 transcripts that were annotated in the GO database were categorized into 58 functional groups, including three main GO ontologies: biological processes, cellular components, and molecular functions (Fig. 3). Among these functional groups, the terms "biological regulation", "cellular process", "metabolic process", "cell", "single-organism process", "cell part", "binding", and "catalytic activity" were dominant.
COGs were delineated by comparing protein sequences encoded in the complete genome, representing major phylogenetic lineages 18 . In this study, COG classification was used to further evaluate the completeness of the transcriptome library and the effectiveness of annotation methods. As a result, a total of 15,959 unigenes were mapped into 25 different COG categories. Of these categories, the largest COG group was the R category, representing "general function prediction only" (6,435 unigenes, 41.26%); followed by the J category, representing "translation, ribosomal structure and biogenesis" (3,347 unigenes, 21.46%); the K category, representing "transcription" (3,044 unigenes, 19.52%); and the L category, representing "replication, recombination and repair" (2,802 unigenes, 17.97%) (Fig. 4).
In addition, KEGG is a bioinformatics resource for linking genomes to life and the environment 19 . The KEGG PATHWAY database records networks of molecular interactions in the cells and variants specific to particular organisms 20 . The genes from the merged groups (GN and GW) were categorized using the KEGG database to obtain more information to predict unigene function 21 . As a result, a total of 25,290 unigenes were classified into 257 KEGG pathways. Among these KEGG pathways, the top 50 statistically significant KEGG classifications are shown in Table 3. Some important innate immunity-related pathways were predicted in this KEGG database, including Vibrio cholerae infection (1092 sequences, 4.32%), focal adhesion (910 sequences, 3.6%), Epstein-Barr virus infection (860 sequences, 3.40%), lysosome (610 sequences, 2.41%), HTLV-I infection (596 sequences, 2.36%), Herpes simplex infection (593 sequences, 2.34%), salmonella infection (576 sequences, 2.28%), MAPK signaling pathway (542 sequences, 2.14%), adherens junction (467 sequences, 1.85%), and so on (Table 3). It is worth noting that the insulin signaling pathway, the Wnt signaling pathway, the mRNA surveillance pathway, endocytosis, phagosome, ECM-receptor interaction, bacterial invasion of epithelial cells, Fc gamma R-mediated phagocytosis, and tight junction were present in the top 50 statistically significant KEGG classification. Perhaps these signaling pathways and related genes have an important effect on further understanding the antiviral mechanisms of the innate immune system.

Differentially expressed gene analysis in crayfish intestine after WSSV infection. Previous
sequence analysis and annotation for all of the unigenes in the merged group (GN and GW) provided some valuable information to analyze the crayfish intestine transcriptome. However, the variation in the gene expression level of crayfish intestine after WSSV infection was expected. In this study, FDR ≤ 0.001 and an absolute value of log2Ratio ≥ 1 were used as the filtering thresholds to identify up-regulated or down-regulated genes between normal crayfish and WSSV-challenged crayfish. As shown in Fig. 5, a total of 7,000 DEGs were screened out after a comparative analysis between the GN and GW samples. Among these genes, 5,976 were identified as differentially up-regulated and 1,024 as differently down-regulated by more than two fold. Among these 7,000 DEGs, 6,821 genes existed in both the GN and GW samples at the same time, including 5,798 differentially up-regulated genes and 1,024 differently down-regulated genes. Moreover, 178 genes were only found in the GW sample after WSSV challenge. Compared to the abovementioned signaling pathways, these 178 genes have an important influence on further studies of the antiviral immune mechanisms.
To determine the biological function of DEGs between GN and GW, GO classification and KEGG pathway analysis were carried out 22 . GO classification analysis was performed on annotated transcripts using Blast2GO. As shown in Fig. 6, a total of 1,270 DEGs were screened out after a comparison between the GN and GW samples. The results showed that 1,270 DEGs that were annotated in the GO database were categorized into 52 functional groups, including the three main GO ontologies: biological processes, cellular components, and molecular functions. Among these DEGs, a large number were dominant in nine terms, including biological regulation, cellular process, metabolic process, single-organism process, cell, cell part, organelle, binding, and catalytic activity.
Then, all of the DEGs were mapped in the KEGG database to search for genes involved in the innate immune response or signaling pathways. A total of 7,000 DEGs were assigned to 250 KEGG pathways. KEGG pathway analysis showed that 36 pathways were obviously changed (P-value < 0.05) in the GW sample compared with the GN sample. Among these 36 pathways, 28 were significantly changed (P-value < 0.01), and some were related to the innate immunity response, including the insulin signaling pathway, the Wnt signaling pathway, ECM-receptor interaction, the JAK-STAT signaling pathway, cell adhesion molecules, the mRNA surveillance pathway, cytokine-cytokine receptor interaction, lysosome, adherens junction, and the Notch signaling pathway (Table 4). Because the antiviral immune mechanism of crayfish is not clear, the discovery of these above signaling pathways for DEGs will help to identify the innate immune mechanisms.
In-depth analysis of DEGs involved in signaling pathways related to innate immunity. Based on the KEGG pathway analysis of DEGs between GN and GW, some classical signaling pathways that were related to the innate immune system were screened out, for example, the JAK-STAT signaling pathway, insulin signaling pathway, Wnt signaling pathway, mRNA surveillance pathway, and Notch signaling pathway. To date, research results regarding the antiviral immune mechanisms of crustaceans have revealed that the JAK-STAT signaling pathway is involved in the antiviral innate immune response of shrimp 23 . However, the roles of the other four signaling pathways in the crustacean antiviral immune response have not been reported.
The JAK-STAT signaling pathway has also been implicated in the insect antiviral immune defense response, which includes three main cellular components: receptor Domeless, Janus Kinase (JAK) Hopscotch, and STAT transcription factor 24 . Moreover, the transcription of STAT in shrimp was obviously up-regulated after WSSV infection 25 , indicating that the JAK-STAT pathway might play a very important role in shrimp antivirus immunity responses. According to the transcriptome sequencing results in the present study, some unigenes were annotated in the JAK-STAT signaling pathway, and their expression levels obviously varied after WSSV infection. These results suggest that this pathway plays an important role in the crayfish antiviral innate immune response. In the present study, 21 genes were significantly differentially expressed in the JAK-STAT signaling pathway, including 19 significantly up-regulated genes and 2 significantly down-regulated genes (Fig. 7). The protein identification and concrete expression profile analysis of these 21 genes is shown in Table 5. Some important molecules involved in the classical JAK-STAT signaling pathway are included, for example, STAT, suppressor of cytokine signaling-2 like protein (SOCS), apoptosis regulator Bcl-XL, Myc protein, Ras GTP exchange factor, and Src  kinase-associated phosphoprotein 2 (SKAP). Based on the transcriptome sequencing results, the gene expression levels of these abovementioned molecules were obviously up-regulated after WSSV infection (Table 5).
To further ascertain the role of the JAK-STAT signaling pathway in crayfish antiviral immune responses, five key DEGs involved in the JAK-STAT signaling pathway were selected for qRT-PCR to analyze their expression profiles in crayfish intestine after WSSV infection. The protein identities of the five DEGs were STAT, Ras GTP exchange factor, Ras GTP exchange factor, apoptosis regulator Bcl-XL, and suppressor of cytokine signaling-2 like protein. The qRT-PCR results showed that these five DEGs were obviously up-regulated at 36 h post WSSV infection in crayfish intestine (Fig. 8). These results could provide new insight into crayfish antiviral immunity. To clarify the functions of this pathway, other components need to be identified, and the interaction among these components needs to be explored as soon as possible.

Validation of transcriptome data by qRT-PCR.
We selected thirteen genes that are related to the innate immunity response to evaluate their differential expression levels between GN and GW samples using qRT-PCR 26 . For these candidate genes, the qRT-PCR expression profile patterns were consistent with the RNA-Seq data ( Table 6). There were similar trends of gene up/down-regulation between the qRT-PCR and data. The results illustrated that the RNA-Seq data were reliable.

Figure 5. Comparative analysis of gene expression levels for two transcript libraries between the normal crayfish intestine (GN) and WSSV-challenged crayfish intestine (GW) samples.
Red dots represent transcripts that were significantly up-regulated in GW, and green dots indicate that those transcripts were significantly down-regulated. The parameters "FDR ≤ 0.001" and "|log2 Ratio| ≥ 1" were used as the thresholds to judge the significance of gene expression differences.

Conclusions
In this study, the de novo-assembled transcriptomes of crayfish intestines were analyzed, and a large amount of sequence information was obtained. The expression profiles of DEGs between the normal crayfish intestine transcriptome and that of WSSV-challenged crayfish was studied. The aim of this deep analysis of DEG functional annotation, orthologous protein clustering, and annotation of signaling pathways related to the immune system was to determine the underlying mechanisms involved in the anti-WSSV immune response in crayfish. Based on the transcriptome sequencing results in the present study, many genes and pathways related to innate immunity in crayfish intestine were regulated after WSSV challenge. It is worth noting that 7,000 DEGs were screened out after a comparative analysis between the GN and GW samples. These DEGs were mapped into 250 KEGG pathways. Among these pathways, 36 were obviously changed (P-values < 0.05) and 28 pathways were extremely significantly changed (P-values < 0.01).
To further identify the signaling pathways that were related to the crayfish antiviral immune response, five key DEGs involved in the JAK-STAT signaling pathway were selected for qRT-PCR. The results showed that all five of these DEGs were obviously up-regulated at 36 h post WSSV infection in crayfish intestine. Taken together, these results provide new insight into the crayfish antiviral immunity mechanism. In addition, these results could also provide an important theoretical basis for solving viral disease problems in crayfish breeding.

Materials and Methods
Preparation of crayfish tissues and immune challenge. P were originally cultured in water tanks at 26-28 °C for at least 5 days and fed twice daily with artificial food throughout the experiment 27 . For WSSV infection, WSSV (3.2 × 10 7 particles per crayfish) was injected into the abdominal segment of each crayfish 1,28,29 . Then, 36 h after challenge, the intestines were collected from no fewer than ten WSSV-challenged crayfish. These samples constituted the WSSV group (GW). The intestines were also collected from at least ten normal crayfish and frozen immediately in liquid nitrogen. These samples constituted the normal group (GN). Then, these two sets of samples were temporarily stored at − 80 °C for total RNA extraction 27 .

RNA isolation and Illumina sequencing.
The two sets (GN and GW) of intestine tissue samples that had been frozen in liquid nitrogen were delivered to the Beijing Genomics Institute-Shenzhen (BGI, Shenzhen, China) for total RNA extraction. In brief, the total RNA from the crayfish intestines was extracted with TRIzol    Transcriptome de novo assembly and analysis. Transcriptome de novo assembly for the two intestine sample (GN and GW) sets was carried out by the RNA-Seq de novo assembly program Trinity 30 . In brief, the raw reads that were generated by the Illumina HiSeq ™ 2000 sequencer were originally trimmed by removing the adapter sequences. After the low-quality reads with quality scores of less than 20 and short reads with lengths of less than 10 bp were removed, high-quality clean reads were obtained to perform transcriptome de novo assembly using Trinity software with the default parameters. Generally, there were three steps, including Inchworm, Chrysalis and Butterfly 31 . In the first step, high-quality clean reads were processed by Inchworm to form longer fragments, which were called contigs. Then, these contigs were connected by Chrysalis to obtain unigenes that could not be extended on either end. These unigenes resulted in de Bruijn graphs. Finally, the de Bruijn graphs were processed by Butterfly to obtain transcripts 32 .   Transcriptome annotation and gene ontology analysis. After transcriptome de novo assembly, the transcripts were used for annotation, including protein functional annotation, COG functional annotation, GO functional annotation, and pathway annotation. These processes are based on sequence similarity with known genes. In detail, the assembled contigs were annotated with sequences available in the NCBI database using the BLASTx and BLASTn algorithms 33 . The unigenes were aligned by a BLASTx search against the protein databases of NCBI, including Nr, Swiss-Prot, KEGG, and COG 34 . Meanwhile, none of the BLASTx hits were aligned by a BLASTn search against the NCBI Nt database. All of the above alignments were executed to establish the homology of sequences with known genes (with a cutoff E-value ≤ 10 −5 ) 35 . Then, the best alignment results were used to determine the sequence direction and protein-coding-region prediction (CDS) of the unigenes. Functional annotation was executed with GO terms (www.geneontology.org) that were analyzed using Blast2GO software (http://www.blast2go.com/b2ghome) 36 . Based on the KEGG database, the complex biological behavior of the genes was analyzed through pathway annotation.

Identification of differentially expressed genes.
To acquire the expression profiles for transcripts in crayfish intestines, cleaned reads were first mapped to all of the transcripts using Bowtie software 37 . Then, DEGs were obtained based on the number of fragments per kilobase of exon per million fragments mapped (FPKM) of the genes, followed by a False Discovery Rate (FDR) control to correct for the P-value 38 . DEGs were identified using EDGER software (empirical analysis of digital gene expression data in R) 39 . For this analysis, the filtering threshold was set as an FDR control of 0.5. Lastly, FDR ≤ 0.001 and the absolute value of log2Ratio ≥ 1 were used as the filtering thresholds to determine the significance of the differentially expressed genes 40 . The justification for using |log2Ratio| ≥ 1 as the filtering threshold was to reduce the statistical workload while obtaining more meaningful differentially expressed genes. Using this method, the differentially expressed genes were identified between GW and GN through a comparative analysis of the above data.
Quantitative real-time PCR validation. Quantitative real-time PCR (qRT-PCR) methods were used to determine the RNA levels for fifteen selected genes that were related to the innate immune response 41 . For qRT-PCR analysis, cDNA templates from the two intestine sample (GN and GW) sets were diluted 20-fold in nuclease-free water and were used as templates for PCR. Gene-specific primer sequences were carefully designed using Primer Premier 6 software based on the sequence of each gene that was identified from the transcriptome library 42 . The specific primers, namely Pc-18 S RNA-qRT-F (5′-tct tct tag agg gat tag cgg-3′) and Pc-18 S RNA-qRT-R (5′-aag ggg att gaa cgg gtt a-3′), were used to amplify the 18S RNA gene as the inner control. qRT-PCR was performed following the manufacturer's instructions for SYBR Premix Ex Taq (Takara, Japan) using a real-time thermal cycler (Bio-Rad, USA) in a total volume of 10 μ l containing 5 μ l of 2× Premix Ex Taq, 1 μ l of the 1:20 diluted cDNA, and 2 μ l (1 μ M) each of the forward and reverse primers. The amplification procedure comprised an initial denaturation step at 95 °C for 3 min, followed by 40 cycles of 95 °C for 15 s and 59 °C for 40 s and melting from 65 °C to 95 °C. Three parallel experiments were performed to improve the integrity of the work 43 . Furthermore, the differentially expressed levels of the target genes (between the GN and GW samples) were calculated by the 2 −ΔΔCT analysis method as described in a previous study 44 . The obtained data were subjected to statistical analysis, followed by an unpaired sample t-test. A significant difference was accepted at a P-value < 0.05. An extremely significant difference was accepted at P < 0.01.