Immune-Related Functional Differential Gene Expression in Koi Carp (Cyprinus carpio) after Challenge with Aeromonas sobria

In order to understand the molecular basis underlying the host immune response of koi carp (Cyprinus carpio), Illumina HiSeqTM 2000 is used to analyze the muscle and spleen transcriptome of koi carp infected with Aeromonas sobria (A. sobria). De novo assembly of paired-end reads yielded 69,480 unigenes, of which the total length, average length, N50, and GC content are 70,120,028 bp, 1037 bp, 1793 bp, and 45.77%, respectively. Annotation is performed by comparison against various databases, yielding 42,229 (non-redundant protein sequence (NR): 60.78%), 59,255 (non-redundant nucleotide (NT): 85.28%), 35,900 (Swiss-Prot: 51.67%), 11,772 (clusters of orthologous groups (COG): 16.94%), 33,057 (Kyoto Encyclopedia of Genes and Genomes (KEGG): 47.58%), 18,764 (Gene Ontology (GO): 27.01%), and 32,085 (Interpro: 46.18%) unigenes. Comparative analysis of the expression profiles between bacterial challenge fish and control fish identifies 7749 differentially expressed genes (DEGs) from the muscle and 7846 DEGs from the spleen. These DEGs are further categorized with KEGG. Enrichment analysis of the DEGs and unigenes reveals major immune-related functions, including up-regulation of genes related with Toll-like receptor signaling, complement and coagulation cascades, and antigen processing and presentation. The results from RNA-Seq data are also validated and confirmed the consistency of the expression levels of seven immune-related genes after 24 h post infection with qPCR. Microsatellites (11,534), including di-to hexa nucleotide repeat motifs, are also identified. Altogether, this work provides valuable insights into the underlying immune mechanisms elicited during bacterial infection in koi carp that may aid in the future development of disease control measures in protection against A. sobria.


Introduction
Aeromonas (Aeromonadaceae) species are ubiquitous and can cause infections not only in humans but also in fish. Aeromonas are isolated from various sources, such as fresh, estuarine, or surface waters,

Number of Differentially Expressed Genes after Aeromonas Sobria Challenge
Comparison of gene expression levels between the fish, subjected to bacterial challenge and control fish, identified a total of 7749 differentially expressed genes (DEGs) in the muscle and 7846 in the spleen (p < 0.05). This included 6300 up-regulated and 1449 down-regulated genes in the muscle, and 5111 up-regulated and 2735 down-regulated genes in the spleen, revealing a total of 7192 and 7280 ( Figure 1). The DEGs in the muscle and spleen were mainly annotated into the following categories: 'biological process', 'cellular component', and 'molecular function' (Figures 2 and 3). The most annotated unigenes belonged to the following categories: cellular process, single-organism process and metabolic process (from the 'biological process' category); cell, cell part and organelle (from the 'cellular component' category); and binding, catalytic activity and molecular transducer activity (from the 'molecular function' category).
(from the 'cellular component' category); and binding, catalytic activity and molecular transducer activity (from the 'molecular function' category).    (from the 'cellular component' category); and binding, catalytic activity and molecular transducer activity (from the 'molecular function' category).    Overall, from DEGs in muscle and spleen, KEGG analysis annotated the genes into 286 pathways, which were classified into 6 main categories (Figures 2 and 3), namely cellular processes, environmental information processing organismal systems, genetic information processing metabolism, human diseases, metabolism, and organismal systems.
The immune system showed the highest number of DEGs in the muscle (817 genes); these were divided into 16 sub-categories (Supplementary Materials, Table S3a,b), including leukocyte transendothelial migration (203 genes), platelet activation (184 genes), hematopoietic cell lineage (127 genes), T cell receptor signaling pathway (147 genes), natural killer cell-mediated cytotoxicity (119 genes), NOD-like receptor signaling pathway (98 genes), Fc gamma R-mediated phagocytosis (145 genes), B cell receptor signaling pathway (107 genes), intestinal immune network for IgA production (55 genes), antigen processing and presentation (75 genes), Toll-like receptor signaling pathway ( Figure          DEGs with an absolute value of fold change >1 were selected from the immune-related category from the muscle and spleen and are presented in Table 1. Although most of the selected genes could be found in other categories, the genes were related to the complement system, antigen processing and presentation, and the Toll-like receptor signaling pathway.  DEGs with an absolute value of fold change >1 were selected from the immune-related category from the muscle and spleen and are presented in Table 1. Although most of the selected genes could be found in other categories, the genes were related to the complement system, antigen processing and presentation, and the Toll-like receptor signaling pathway. DEGs with an absolute value of fold change >1 were selected from the immune-related category from the muscle and spleen and are presented in Table 1. Although most of the selected genes could be found in other categories, the genes were related to the complement system, antigen processing and presentation, and the Toll-like receptor signaling pathway.

RT-qPCR Analysis of Immune Related Genes Following Aeromonas sorbia Infection
To validate the DEGs identified by RNA-Seq analysis, we performed RT-qPCR. Seven genes (C3, IL1β, IL8, MyD88, NF-κB, TLR5, TNFα) associated with the complement system, antigen processing, and toll-like receptors were detected in the spleen and muscle by RT-qPCR (Table 4). In the spleen of the infected fish, at 1 dpi, the expression levels of IL1β, IL8, and MyD88 were significantly up-regulated. In the muscle at 1 dpi, the expression levels of C3, IL1β, and IL-8 showed an upward trend, but these was not significant.

Discussion
This study used the spleen and muscle of koi carp at 24 h after infection with A. sobria as experimental samples to determine immune-related genes and signaling pathways activated during the early stage of infection. As expected, the results showed that many immune-related genes in the koi carp were up-regulated significantly after A. sobria infection. The most significantly up-regulated genes associated with immunity were the pro-inflammatory cytokine-related and the signal transduction related genes, such as IL-1β, TNF receptor, CXC chemokine, TGF-β, NF-κB, and some other immune-related genes; pathogen recognition related genes were also significantly up-regulated.
After assembly, 69,480 unigenes were generated, with an average length of 578 bp, which was longer than those achieved in previous studies, using the Roche GS FLX 454 system with MIRA assembler (a length range of 118-2065 bp and an average length of 495 bp) [32] or Illumina/Hiseq-2000 with the assembling program-SOAP (a length range of 200-5245 bp and an average length of 412 bp) [33]. This difference in sequence quality is possibly owing to the difference in sampling tissue and the different de novo assemblers. Trinity assemblies, which were used in this study, have a consistently better performance than the other tools used in transcriptome assembling, even in the absence of a reference genome [34,35]. In contrast to Trinity, the SOAP or MIRA assemblies adopted in previous reports [32,33] were more fragmented under high values of

Discussion
This study used the spleen and muscle of koi carp at 24 h after infection with A. sobria as experimental samples to determine immune-related genes and signaling pathways activated during the early stage of infection. As expected, the results showed that many immune-related genes in the koi carp were up-regulated significantly after A. sobria infection. The most significantly up-regulated genes associated with immunity were the pro-inflammatory cytokine-related and the signal transduction related genes, such as IL-1β, TNF receptor, CXC chemokine, TGF-β, NF-κB, and some other immune-related genes; pathogen recognition related genes were also significantly up-regulated.
After assembly, 69,480 unigenes were generated, with an average length of 578 bp, which was longer than those achieved in previous studies, using the Roche GS FLX 454 system with MIRA assembler (a length range of 118-2065 bp and an average length of 495 bp) [32] or Illumina/Hiseq-2000 with the assembling program-SOAP (a length range of 200-5245 bp and an average length of 412 bp) [33]. This difference in sequence quality is possibly owing to the difference in sampling tissue and the different de novo assemblers. Trinity assemblies, which were used in this study, have a consistently better performance than the other tools used in transcriptome assembling, even in the absence of a reference genome [34,35]. In contrast to Trinity, the SOAP or MIRA assemblies adopted in previous reports [32,33] were more fragmented under high values of sequencing errors and polymorphism levels [34]. Additionally, compared with the merged 69,480 unigenes in this transcriptome, previous studies provided relatively smaller gene sets (29,682 unigenes from Roche 454 system and 2139 unigenes from a SMART cDNA library) [36], and this result further emphasizes that Illumina/Hiseq-2000 RNA-Seq is a more ideal method for transcriptome analysis and it has high efficiency and massive data output.
In this study, only 60,593 unigenes were annotated with at least one database and the same problem has also been reported for the transcriptomes of other groups of marine organisms [37]. The possible reason for less annotation may be due to the availability of limited genome sequence database and research on aquaculture fish species [38][39][40]. Nonetheless, the functional annotations of unigenes according to GO, COG, and KEGG databases provided ample numbers of candidate genes and valuable information about biological features of koi carp challenged with A. sobria in this study. For example, as per the KEGG analysis, 33,057 sequences were assigned to 244 KEGG pathways, and among them, genetic information processing accounted for the largest number of pathways related to pathogen infection (Figures 2 and 3).
In the present study, identification of SSR might be useful in genetic, evolutionary, and breeding studies [41] in koi carp studies. These data on SRR will be useful in future studies on gene expression, which can be manipulated by SRR variations in the 5 UTR regions as they affect the transcription and translation. As it is evident that SRR are ubiquitous in transcriptomes and specific, this unigene could be used for molecular marker development, comparative genetic mapping and genotyping.

Animal Maintenance
Prior to the conduct of experiments, all the fish were kept in a recirculatory system for 2 weeks and allowed them to acclimatize to the laboratory conditions. Throughout the experiment, fish were handled with 2-phenoxyethanol as an anesthetic. Approval for animal studies was obtained from the Center for Research Animal Care and Use Committee of the National Pingtung University of Science and Technology under protocol no #101-027 dated 19 March 2012.

Isolation, Cultivation, and Challenge with Aeromonas Sobria
The bacterium A. sobria was isolated from diseased koi carp with severe skin ulcer. The species was identified by API 20NE and 16S rDNA sequencing, cultured in brain heart infusion (BHI) broth for 18 h at 25 • C and enumerated prior to the challenge test.
A total of 15 fish (body weight 325 ± 23 g) were anaesthetized and used for intraperitoneal injection with 1 × 10 7 cfu per fish. Individual fish received 1 × 10 7 cfu in 200 µL phosphate buffered saline (PBS, pH 7.2). The other tank with 15 fish was injected with 200 µL of PBS (pH 7.2) only and used as a control. Three fish each from the challenge (treatment) and control groups (n = 3), respectively, were examined at post 24 h infection. Muscle and spleen tissues were dissected hygienically and total RNA was isolated.

Total RNA Isolation
TRIzol ® reagent (Invitrogen Corp., Carlsbad, CA, USA) was used to isolate total RNA from the tissues according to the manufacturer's instructions. RNA integrity was assessed using the RNA Nano 6000 Assay Kit on the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

cDNA Library Preparation and Sequencing
Genomics Bioscience Technology Co. Ltd. (Taipei, Taiwan) synthesized cDNA using 40 µg total RNA along with poly-T oligo-attached magnetic beads. First-and second-strand cDNA was synthesized using random oligonucleotides and SuperScript II reverse transcriptase. Illumina HiSeq™ 2000 (Illumina, Inc., San Diego, CA, USA) platform was used to sequence the RNA-Seq library as paired-end reads to 100 bp.

Assembly of De Novo Transcriptome
Adaptors and unknown bases (N) and low-quality reads were filtered using internal software and stored in FASTQ format [42]. Using the default parameter in Trinity (https://github.com/ trinityrnaseq/trinityrnaseq/wiki) assembly was performed for paired reads and clustered to unigenes via TIGR gene indices.

Differentially Expressed Genes and Enrichment Analysis
Bowtie2 software (http://bowtie-bio.sourceforge.net) [49] was used to determine the expression form treatment and control library. In order to obtain the DEGs in muscle and spleen tissues between the control and the infected groups, fragments per kilobase of transcripts per million fragments mapped (FPKM) values were analyzed further using the RESM [50] and the false discovery rate (FDR) was used when it is <0.05.

Real-Time Reverse Transcription Polymerase Chain Reaction
Once the infected fish were prepared, we performed another set of experiments for the validation of RNA-Seq data. DNase I-treated total RNA from the spleen and muscle was subjected to cDNA synthesis using iScript™ cDNA synthesis kits (Bio-Rad, Hercules, CA, USA). Reverse transcriptase real-time PCR (RT-qPCR) was performed using iQ™ SYBR ® Green Supermix (Bio-Rad). The list of primer sequences is shown in Table 5. Gene expression levels were normalized to that of EF1α. To enable comparisons between the two groups, statistical analysis was performed using Student's t-test. Values of p < 0.05 were considered statistically significant. Table 5. Primer name and sequence used in the present study.

Conclusions
This is the first study to provide information on host defense gene activities based on differential transcriptomic profiling in koi carp against A. sobria. The large number of differential expression of immune-related genes from innate immunity, antigen processing and complement cascade in this study indicated the activation of koi carp host immune response towards A. sorbria infection and its effect. However, further studies should be directed towards understanding of these different expressed immune genes as potential functional markers in koi carp against A. sobria infection. Understanding the protein levels of these functional markers in the infected koi carp would be of vital importance, as such approach may lead to a profound understanding of the regulation and responses that occur during the infection, and thus lead to the development of very specific and efficient novel vaccines and chemotherapies.