Comparative transcriptome analysis reveals three potential antiviral signaling pathways in lymph organ tissue of the red swamp crayfish, Procambarus clarkii

Ethics statement

The experimental procedures complied with the current applicable laws of China where they were performed. No specific permits were required for the research undertaken. Crayfish individuals were maintained under appropriate laboratory conditions to guarantee their welfare and responsiveness.

Preparation of crayfish tissues and WSSV challenge

P. clarkii (weighing approximately 15-20 g) were bought from an aquaculture commercial market in Hangzhou, Zhejiang Province, China. The crayfish were initially cultured in water tanks at 26°-28°C for 10 days, and were fed twice a day with artificial food throughout the experiment (Li et al., 2012b). WSSV (3.2 x 10⁷ copies per crayfish) was injected into the abdominal segment of each crayfish; 36 h later, the lymph organ was removed from at least 10 WSSV-challenged crayfish [Group WSSV-challenged (GW)]. The lymph organs were also removed from at least 10 normal crayfish and frozen immediately in liquid nitrogen [Group normal (GN)]. The two sets of samples were temporarily stored at -80°C for total RNA extraction (Du and Jin, 2015).

RNA isolation and Illumina sequencing

The samples were delivered to the Beijing Genomics Institute-Shenzhen (Shenzhen, China) for total RNA extraction. According to the company’s sequencing report, total RNA was extracted using TRIzol reagent (Invitrogen, USA) in accordance with the manufacturer protocol. The quality of RNA samples treated with DNase I (Invitrogen) was examined for the subsequent procedures, including mRNA purification, cDNA library construction, and transcriptome sequencing. In brief, about 5 μg DNase-treated total RNA was used to construct a cDNA library following the protocol of the Illumina TruSeq RNA Sample Preparation Kit (Illumina, USA). After some necessary quantifications and qualifications, the library was sequenced using an Illumina HiSeq^TM 2000 instrument with 100-bp paired-end reads for GN and GW.

De novo assembly and transcriptome analysis

De novo transcriptome assembly of the two sets of samples (GN and GW) was conducted using the RNA-Seq de novo assembly program, Trinity (Dedeine et al., 2015), with the default parameters set (Ali et al., 2015). In brief, raw reads generated by the Illumina HiSeq^TM 2000 sequencer were trimmed by removing the adapter sequences. After low-quality reads (quality scores lower than 20) and short-length reads (below 10 bp) were removed, high-quality clean reads were obtained. Three processes were then performed using Inchworm, Chrysalis, and Butterfly (He et al., 2015). Firstly, high-quality clean reads were processed by Inchworm to form longer fragments called contigs. The contigs were connected by Chrysalis to obtain unigenes that could not be extended at either end; the unigenes resulted in de Bruijn graphs. Finally, the de Bruijn graphs were processed using Butterfly in order to obtain transcripts (Li et al., 2013).

Transcriptome annotation and Gene Ontology (GO) analysis

After the de novo transcriptome assembly was completed, the transcripts were used to conduct annotation (Mousavi et al., 2014), including protein functional annotation, Cluster of Orthologous Groups (COG) functional annotation, GO functional annotation (www.geneontology.org), and pathway annotation. All of these were based on sequence similarity with known genes. The assembled contigs were annotated with sequences available in the National Center for Biotechnology Information (NCBI) database using the BLASTx and BLASTn algorithms (Leng et al., 2015). The unigenes were aligned by a BLASTx search against NCBI protein databases, including the non-redundant sequence database (Nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and COG (Li et al., 2015). None of the BLASTx hits were aligned by a BLASTn search against the NCBI non-redundant nucleotide database (Nt). The alignments were conducted in order to establish the homology of sequences with known genes with a cutoff E-value ≤ 10^-5 (Li et al., 2012a). The best alignment results were used to predict the sequence direction and the coding DNA sequence of the unigenes. Functional annotation was conducted using GO terms that were analyzed using the Blast2go software (http://www.blast2go.com/b2ghome). The complex biological behavior of the genes was analyzed by pathway annotation based on the KEGG database.

DEG identification

To ascertain transcript expression levels in the crayfish lymph organ, clean reads were first mapped to all of the transcripts using the Bowtie software (Langmead and Salzberg, 2012). DEGs were then obtained on the basis of fragments per kilobase of exon per million fragments mapped (FPKM) of the genes, followed by a false discovery rate (FDR) control to correct the P values (Mortazavi et al., 2008). DEGs were identified using the EDGER software. For this analysis, the filtering threshold was set as an FDR control, 0.5. FDR ≤ 0.001 and the absolute value of log2Ratio ≥ 1 were used as the filtering threshold to determine DEG significance (Banani et al., 2016). DEGs were identified between GW and GN by a comparative analysis of the data.

Quantitative real-time polymerase chain reaction (PCR) validation

Quantitative real-time PCR (qRT-PCR) methods were used to determine the RNA levels of 10 randomly selected genes (Wang et al., 2016). cDNA templates from the two sets of samples (GN and GW) were diluted 30-fold in nuclease-free water and used as templates for the PCR. Gene-specific primer sequences were carefully designed using the Primer Premier 6 software, based on each gene sequence identified in the transcriptome library. The specific primers 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 18 S RNA as an inner control. qRT-PCR was performed using SYBR^® Premix Ex Taq (TaKaRa, Japan) and a real-time thermal cycler (CFX96 touch, Bio-Rad, USA) following the manufacturer instructions, in a total volume of 10 μL that contained 5 μL 2X Premix Ex Taq, 1 μL 1:30 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 at 95°C for 15 s, 60°C for 40 s, and melting from 65° to 95°C. Three parallel experiments were performed to improve the integrity of the study. Differentially expressed levels of the target genes (between GW and GN) were calculated using the 2^-ΔΔCt method (Livak and Schmittgen, 2001). The data obtained were subjected to an unpaired-samples Student t-test. Significance was set at P < 0.05.

Illumina sequencing of the crayfish lymph organ

After cleaning and quality testing the GN samples, 44,264,442 clean reads were screened from 46,182,460 raw reads, which corresponded to 4,426,444,200 total nucleotides (nt). The Q20 (percentage of bases the quality of which was greater than 20 in clean reads), N (percentage of uncertain bases after filtering), and GC percentages were 97.92, 0.00, and 42.99%, respectively. For the GW samples, 50,263,784 clean reads were screened from 52,010,554 raw reads, which corresponded to 5,026,378,400 total nt. The Q20, N, and GC percentages were 97.86, 0.00, and 43.68%, respectively. All of these sequences were used for further analysis.

De novo assembly of the transcriptome

For the GN samples, 44,264,442 high-quality clean reads were randomly assembled to produce 105,034 contigs with an N50 of 546 bp. The contigs were further assembled and clustered into 56,922 unigenes with a mean length of 564 nt, which were composed of 4886 distinct clusters and 52,036 distinct singletons. The N50 of the unigenes was 1032 bp. For the GW samples, 50,263,784 high-quality clean reads were randomly assembled to produce 137,886 contigs with an N50 of 611 bp. The contigs were further assembled and clustered into 74,732 unigenes with a mean length of 606 nt, which were composed of 7646 distinct clusters and 67,086 distinct singletons. The N50 of the unigenes was 1197 bp.

In the GN samples, the unigenes were mainly between 200 and 3000 nt long. The percentage of unigenes that were between 201 and 1000 bp long was 86.9% (49,456); between 1001 and 2000 bp long it was 8.0% (4547); between 2001 and 3000 bp long it was 2.7% (1538); and for over 3000 bp long it was 2.4% (1381). In the GW samples, the percentage of unigenes that were between 201 and 1000 bp long was 85.5% (63,916); between 1001 and 2000 bp long it was 8.5% (6337); between 2001 and 3000 bp long it was 3.1% (2303); and for over 3000 bp long it was 2.9% (2176).

Functional annotation of predicted proteins

Sequence annotation was conducted based on the unigenes from the merged group. Subsequently, the putative functions of the unigenes were analyzed based on their GO and COG classifications (Young et al., 2010). GO classification is an international, standardized gene function classification system, which provides a dynamically updated controlled vocabulary and an exactly defined conception to describe gene characteristics and their products in any organism (Di Lena et al., 2015). GO classification consists of three main categories: “biological process”, “cellular component”, and “molecular function” (Chen et al., 2015a). A total of 13,848 transcripts annotated in the GO database were categorized into 58 functional groups, including the three main GO categories. Among these functional groups, the “biological regulation”, “cellular process”, “metabolic process”, “cell”, “single-organism process”, “cell part”, “binding”, and “catalytic activity” terms were dominant.

COGs were delineated by comparing protein sequences encoded in complete genomes that represented major phylogenetic lineages (Tatusov et al., 2003). COG classification analysis is an important method of unigene functional annotation and evolutionary research (Lai and Lin, 2013). A total of 37,531 unigenes were mapped to 25 COG categories (Figure 1). The largest COG group was the R category, which represented “general function prediction only” (6435 unigenes), followed by the J category that represented “translation, ribosomal structure, and biogenesis” (3347 unigenes), the K category that represented “transcription” (3044 unigenes), and the L category that represented “replication, recombination, and repair” (2802 unigenes).

Figure 1

Cluster of Orthologous Groups (COG) classification of putative proteins.

KEGG is a bioinformatic resource for linking genomes to life and the environment (Mao et al., 2005). The PATHWAY database record networks of molecular interactions in cells, and variants of them that are specific to particular organisms (Chen et al., 2015b). The genes from the merged groups (GN and GW) were categorized using the KEGG database in order to predict the unigenes’ functions (Feng et al., 2015). A total of 25,290 unigenes were classified into 257 KEGG pathways, the 25 most statistically significant of which are shown in Table 1. Important innate immunity-related pathways were predicted, including Vibrio cholerae infection (1092 sequences), focal adhesion (910 sequences), Epstein-Barr virus infection (860 sequences), HTLV-I infection (596 sequences), Herpes simplex infection (593 sequences), Salmonella infection (576 sequences), lysosome (610 sequences), and the MAPK signaling pathway (542 sequences). The mRNA surveillance pathway, endocytosis, phagosome, and tight junction were amongst the 25 most statistically significant KEGG classifications.

DEG analysis of the crayfish lymph organ after WSSV infection

As shown in Figure 2, 5933 DEGs were screened after a comparative analysis of the GN and GW samples; 4638 were differentially upregulated and 1295 were differentially downregulated by more than two-fold. Of the 5933 DEGs, 5603 were in both the GN and GW groups, including 4393 differentially upregulated and 1210 differentially downregulated genes. In addition, 245 genes were only found in the GW group and 85 were only found in the GN group.

Figure 2

Comparative analysis of gene expression levels for two transcript libraries of normal crayfish lymph organ (GN) and white spot syndrome virus-challenged crayfish lymph organ (GW). Red dots represent transcripts that were significantly upregulated in GW and green dots indicate transcripts that were significantly downregulated. The parameters ‘‘FDR ≤ 0.001’’ and ‘‘log2Ratio ≥ 1’’ were used as the criteria by which to judge the significance of the gene expression difference. DEG, differentially expressed gene; FPKM, fragments per kilobase of exon per million fragments mapped.

In order to confirm the biological functions of the DEGs in the GN and GW groups, GO classification and KEGG pathway analysis were conducted (Gupta et al., 2015). GO analysis of the 5933 DEGs revealed that 47 GO terms were significantly expressed (either up or down; P < 0.05) and 36 were highly significantly expressed (either up or down; P < 0.01). Among the 47 GO terms that had P values lower than 0.05, 31, 7, and 9 belonged to “cellular component”, “molecular function”, and “biological process”, respectively. Among the 36 GO terms that had P values lower than 0.01, 25, 5, and 6 belonged to “cellular component”, “molecular function”, and “biological process”, respectively. Gene clusters with significant expression differentiation were mainly related to “cellular component”.

The KEGG pathway analysis revealed that 35 pathways were significantly different (P < 0.05) between the groups. Of these, 24 were highly significantly different between the groups (P < 0.01). Some were related to crayfish innate immunity, including the JAK-STAT signaling pathway, insulin signaling pathway, mRNA surveillance pathway, Wnt signaling pathway, extracellular matrix-receptor interaction, cell adhesion molecules, peroxisome, lysosome, RNA degradation, and basal transcription factors (Table 2).

In-depth analysis of DEGs involved in antiviral innate immunity signaling pathways

Based on the KEGG pathway analysis of the DEGs between GW and GN, some signaling pathways that were related to the invertebrate innate immune system were screened out, including the JAK-STAT, insulin, and Wnt signaling pathways. All three signaling pathways were significantly different between the groups (P < 0.01). In the JAK-STAT signaling pathway, JAK activation stimulates cell proliferation, differentiation, migration, and apoptosis (Harel et al., 2015). Cell proliferation and apoptosis are involved in the antiviral immune response (Du et al., 2013). Eighteen innate immunity-related genes were significantly differentially expressed in the JAK-STAT signaling pathway, 15 of which were significantly upregulated and 3 significantly downregulated (Figure 3). In the insulin signaling pathway, 38 innate immunity-related genes were significantly differentially expressed, including 29 upregulated and 9 downregulated genes (Figure 4). In the Wnt signaling pathway, 19 innate immunity-related genes were significantly differentially expressed, including 15 upregulated and 4 downregulated genes (Figure 5).

Figure 3

Significantly differentially expressed genes identified by the Kyoto Encyclopedia of Genes and Genomes database as being involved in the JAK-STAT signaling pathway. Red boxes indicate significantly increased expression, green boxes indicate significantly decreased expression, and blue boxes indicate unchanged expression.

Figure 4

Significantly differentially expressed genes identified by the Kyoto Encyclopedia of Genes and Genomes database as being involved in the insulin signaling pathway. Red boxes indicate significantly increased expression, green boxes indicate significantly decreased expression, and blue boxes indicate unchanged expression.

Figure 5

Significantly differentially expressed genes identified by the Kyoto Encyclopedia of Genes and Genomes database as being involved in the Wnt signaling pathway. Red boxes indicate significantly increased expression, green boxes indicate significantly decreased expression, and blue boxes indicate unchanged expression.

Validation of transcriptome data by qRT-PCR

We randomly chose 10 genes that were related to the immune system to evaluate their expression levels in the GW and GN samples using qRT-PCR (Gao et al., 2015). The results of the qRT-PCR were consistent with those obtained from the RNA-Seq data (Table 3) and had similar up/downregulation patterns, suggesting that the RNA-Seq data were reliable.