Characterization of the Transcriptome and Gene Expression of Brain Tissue in Sevenband Grouper (Hyporthodus septemfasciatus) in Response to NNV Infection

Grouper is one of the favorite sea food resources in Southeast Asia. However, the outbreaks of the viral nervous necrosis (VNN) disease due to nervous necrosis virus (NNV) infection have caused mass mortality of grouper larvae. Many aqua-farms have suffered substantial financial loss due to the occurrence of VNN. To better understand the infection mechanism of NNV, we performed the transcriptome analysis of sevenband grouper brain tissue, the main target of NNV infection. After artificial NNV challenge, transcriptome of brain tissues of sevenband grouper was subjected to next generation sequencing (NGS) using an Illumina Hi-seq 2500 system. Both mRNAs from pooled samples of mock and NNV-infected sevenband grouper brains were sequenced. Clean reads of mock and NNV-infected samples were de novo assembled and obtained 104,348 unigenes. In addition, 628 differentially expressed genes (DEGs) in response to NNV infection were identified. This result could provide critical information not only for the identification of genes involved in NNV infection, but for the understanding of the response of sevenband groupers to NNV infection.


Introduction
Grouper is one of the highest valued marine fish and has become an important species in the aquaculture industry of various Asian countries. In Korea, sevenband grouper (Hyporthodus septemfasciatus) is one the favorite grouper fish consumed. Its production rate is increasing. However, viral nervous necrosis (VNN) causes high mortality, especially at the larval and juvenile stage of sevenband groupers during the summer season, which has caused vast economic losses [1]. Viral Nervous Necrosis is a serious disease in the world aquaculture industry [2][3][4]. Firstly, it was reported in bigeye trevally (Caranx sexfasciatus) in the 1980s and since then it has been reported in over twenty species [2][3][4]. The infected fish are usually swimming abnormally and having vacuolization and necrosis of the central nervous system in the brain [3]. In Korea, mass mortalities caused by VNN have been reported from various cultured marine fish such as sevenband grouper (Hyporthodus septemfasciatus), rock bream (Oplegnathus fasciatus), red drum (Sciaenops ocellatus) and olive flounder (Paralichthys olivaceus) since 1990 [5][6][7].
Nervous necrosis virus (NNV), the causative agent of VNN, has non-enveloped icosahedral structure and belongs to the family Nodaviridae (genus Betanodavirus). Its genome contains two single-stranded positive senses RNA: RNA1 (approximately 3.1 kb in length) encodes an RNA-dependent RNA died from day 3 after infection and showed 100% of cumulative mortality after 1 week. The moribund fish at days 3 and 4 were selected for sampling. Brain tissues of three of ten challenged fish from mock and the virus-challenge group were collected and pooled for NGS analysis, respectively.

Next Generation Sequencing of Transcriptome
To obtain high-throughput transcriptome data of sevenband grouper, complementary DNA (cDNA) libraries were prepared for 100 bp paired-end sequencing using a TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocols. They were then paired-end (2 × 100 bp) sequenced using an Illumina HiSeq2500 system (Illumina, San Diego, CA, USA).

Transcriptome Assembly and Functional Annotation
Prior to de novo assembly, paired-end sequences were filtered and cleaned using an NGS QC toolkit [21] to remove low quality reads (Q < 20) and adapter sequences. In addition, bases of both ends less than Q20 of filtered reads were removed additionally. This process is to enhance the quality of reads due to mRNA degradation in both ends of it as time goes on [22]. Only high quality reads were used for de novo assembly performed by Trinity (version 20130225) using default values [23]. To remove the redundant sequences, CD-HIT-EST [24] was used. NCBI Blast (version 2.2.28) was applied for the homology search to predict the function of unigenes. The function of unigenes was predicted by Blastx to search all possible proteins against the NCBI Non-redundant (NR) database (accessed on 17 July 2013). The criterion regarding significance of the similarity was set at Expect-value less than 1 × 10 −5 .

Differentially Expressed Genes Analysis
After obtaining the assembled transcriptome data using Trinity, gene expression level was measured with RNA-Seq by Expectation Maximization (RSEM), a tool for measuring the expression level of transcripts without any information on its reference [25]. The TCC package was used for DEG analysis through the iterative DEGES/DEseq method [26]. Normalization was progressed three times to search meaningful DEGs between comparable samples [27]. The DEGs were identified based on the p-value of less than 0.05.

GO Enrichment of Differentially Expressed Genes
The GO database classifies genes according to the three categories of Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) and provides information on the function of genes. To characterize the identified genes from DEG analysis, a GO based trend test was carried out through the Fisher's exact test. Selected genes with p-values of less than 1 × 10 −5 were regarded as statistically significant.

Data Deposition
All the raw read files were submitted to the sequence reads archive (SRA), NCBI database (accession number-SRR5091816).

The Most Abundantly Expressed Gene in the Transcriptome Profile
To estimate gene expression levels, we calculated the abundance of reads in the transcriptome. The top 20 most highly expressed transcripts are shown in Table 2. Commonly, the most abundant genes in both the mock and NNV-infected groups were ribosomal proteins, such as ribosomal protein (RPS) 15a, RPL39, RPS28, RPS14, RPLS2, RPS27a, RPL21 and RPL32 essential for biological metabolism. Ubiquitin-like protein 4a (UBL4a), C-C motif chemokine 2 (CCL2), lysozyme g (LYG_EPICO) and two novel genes (ID: SGU016297, SGU008676) were highly expressed in the

The Most Abundantly Expressed Gene in the Transcriptome Profile
To estimate gene expression levels, we calculated the abundance of reads in the transcriptome. The top 20 most highly expressed transcripts are shown in Table 2. Commonly, the most abundant genes in both the mock and NNV-infected groups were ribosomal proteins, such as ribosomal protein (RPS) 15a, RPL39, RPS28, RPS14, RPLS2, RPS27a, RPL21 and RPL32 essential for biological metabolism. Ubiquitin-like protein 4a (UBL4a), C-C motif chemokine 2 (CCL2), lysozyme g (LYG_EPICO) and two novel genes (ID: SGU016297, SGU008676) were highly expressed in the NNV-infected group compared to that in the mock group. Of them, Ubiquitin-like protein 4a was the most abundant gene in the NNV-infected group.

Functional Annotations
Putative annotations of these transcripts were performed using BlastX as mentioned in the method section. After gene annotation by using BlastX against the non-redundant (nr) database, the putative functions of 43,280 sequences (41.5%) of 104,348 unigenes sequences were identified.

Immune Relevant DEGs Involved in NNV Infection
A total of 3418 unigenes were differentially expressed based on DEG analysis using the TCC package. A total of 372 genes from the total of 3418 DEGs were annotated (Table S1). Immune relevant DEGs were further manipulated manually (Table 3). In immune relevant genes, a variety of cytokines were intensely up-regulated after NNV infection.
Cathepsins are lysosomal cysteine enzymes with important roles in cellular homeostasis and innate immune response [29]. Among a dozen members of the Cathepsin family, subtypes L, H, K, O, S and Z were up-regulated in the brain of sevenband grouper after NNV infection. Specifically, Cathepsin L was highly expressed in the NNV-infected group showing 8.3 Log FC.
Several lectins were expressed in higher levels in the NNV-infected group compared to the mock group, including C-type lectins (CLEC4M, CLEC10A), galectins (LGALS9, LGALS3), fucolectin (FUCL4), and mannose-binding lectin (MBL). In the case of C-type lectins, its receptor (CD209) was also highly expressed in the infected group (Table 3) indicating that C-type lectin might play specific roles in the response of sevenband grouper to NNV infection.
As expected, a number of antiviral proteins also showed high levels of expression in the NNV-infected group. For example, radical S-adenosyl methionin domain-containing protein 2 (RSAD2), also known as viperin, was highly expressed in the NNV-infected group with 10.40 Log FC.
Mx gene (MX), one of the important downstream effectors of interferon (IFN), was also expressed more in the infected group with 8.64 Log FC. Besides Mx, a lot of IFN-induced proteins were upregulated by NNV infection, including IFN-induced protein 44 (IFI44), IFN-induced protein with tetratricopeptide repeats 5 (IFIT5), IFN-induced very large GTPase 1 (GVINP1), IFN-induced double-stranded RNA-activated protein kinase (EIF2Ak2), and IFN-induced helicase C domain-containing protein 1 (IFIH1). Interestingly, of the various heat shock proteins (HSPs), only HSP30 was significantly upregulated in the NNV-infected group with Log 8.42 FC. NK-Lysin, a known antibacterial protein, was also highly expressed in the NNV-infected group with 8.85 Log FC.

GO Enrichment of Differentially Expressed Genes
GO is a widely used method to classify gene functions and their products in organisms. Therefore, the identified DEGs were subsequently used for GO enrichment analysis. According to GO terms, 2094 (61.3%) of the total of 3418 DEGs were classified into the three categories of molecular function, biological process, and cellular component. "Binding" (1258 genes, 46.3%) was the major subcategory in the molecular function. The largest subcategory found in the biological process category was "cellular process" (1488 genes, 12.3%) while "Cell" (1687 genes, 19.6%) and "cell part" (1687 genes, 19.6%) were the most abundant GO terms in the cellular component category (Figure 2). Because one gene could be categorized into several subcategories, the sum of genes in the subcategories could exceed 100%. GO analysis of the transcriptome revealed nine molecular function subcategories, 62 biological process subcategories, and 12 cellular component subcategories with p value of less than 1 × 10 −5 ) (Table S2).
(61.3%) of the total of 3418 DEGs were classified into the three categories of molecular function, biological process, and cellular component. "Binding" (1258 genes, 46.3%) was the major subcategory in the molecular function. The largest subcategory found in the biological process category was "cellular process" (1488 genes, 12.3%) while "Cell" (1687 genes, 19.6%) and "cell part" (1687 genes, 19.6%) were the most abundant GO terms in the cellular component category (Figure 2). Because one gene could be categorized into several subcategories, the sum of genes in the subcategories could exceed 100%. GO analysis of the transcriptome revealed nine molecular function subcategories, 62 biological process subcategories, and 12 cellular component subcategories with p value of less than 1 × 10 −5 ) (Table S2).

Discussion
NNV infection has caused high mortalities of sevenband groupers in aqua-farms during the summer season, especially at larval and juvenile stages. It has also caused tremendous economic losses [1]. Due to the greater damage to the sevenband grouper industry, investigation on the molecular response of NNV infection is required to understand the outbreak mechanism of disease and develop prevention methods such as vaccines. In this study, we performed a transcriptome analysis of the brain tissue of sevenband grouper infected with NNV compared to mock brain tissue using a RNA-Seq.
Gene annotation by BlastX provides valuable information about the transcripts. In this study, 43,289 unigenes (41.5%) of 104,348 unigenes were annotated. This is similar to the result of orange-spooted grouper (45.8%) [32]. Liu

Discussion
NNV infection has caused high mortalities of sevenband groupers in aqua-farms during the summer season, especially at larval and juvenile stages. It has also caused tremendous economic losses [1]. Due to the greater damage to the sevenband grouper industry, investigation on the molecular response of NNV infection is required to understand the outbreak mechanism of disease and develop prevention methods such as vaccines. In this study, we performed a transcriptome analysis of the brain tissue of sevenband grouper infected with NNV compared to mock brain tissue using a RNA-Seq.
Gene annotation by BlastX provides valuable information about the transcripts. In this study, 43,289 unigenes (41.5%) of 104,348 unigenes were annotated. This is similar to the result of orange-spooted grouper (45.8%) [32]. Liu et al. have addressed the possible reasons of such poor annotation: (1) novel genes; (2) sequencing errors; and (3) artefacts by assembly algorithm [33]. Therefore, more genetic studies are needed to understand the biological functions.
The importance of innate defense mechanisms against viral infection has been extensively reviewed [34][35][36]. In this study, we identified innate immune response relevant genes of sevenband grouper involved in NNV infection. Chemokines are critical components of the immune system. The roles of chemokines and their receptors in viral interactions have been reported in various studies [37]. The chemokines family comprises four subfamilies based on N-terminal cystein-motifs: C, C-C, C-X-C, and C-X3-C subfamilies [38]. In this study, we also detected significant up-regulation of CCL2, CCL34, CCL19, CCL4, CXCL13, CXCL6, CXCL8, and CXCL9 in sevenband grouper brain tissue after NNV infection. Especially, CCL2 was highly over expressed at about 10.66 Log FC. CCL2 is a pro-inflammatory chemokine that is induced during several human acute and chronic viral infections, including human immunodeficiency virus (HIV) [39], hepatitis C virus [40], Epstein-Barr virus [41], respiratory synctitial virus [42], Severe Acute Respiratory Syndrome (SARS) [43], herpes simplex virus-1 [44], and Japanese encephalitis virus [45].
Cathepsins are lysosomal cysteines that play important roles in normal metabolism for the maintenance of cellular homeostasis. Cathepsins are one of the superfamilies involved in the regulation of antigen presentation and degradation as well as immune responses, including apoptosis, inflammation, and regulation of hormone processing [46][47][48]. In addition, Chandran et al. have shown that Cathepsin B and Cathepsin L are involved in Ebola virus infection [49]. They are involved in the entry of reovirus [50]. Recently, Cathepsin L has also been shown to be involved in the entry mediated by the SARS coronavirus spike glycoprotein [51] as well as in the process of fusion glycoprotein of Hendra virus [52]. In this study, Cathepsin L and Cathepsin S were found to be notably expressed after NNV infection. Their functional roles in the interaction between grouper and NNV merit further studies.
Lectins are carbohydrate-binding proteins that are highly specific for sugar moieties. They mediate the attachment and binding of viruses to their targets [53]. Lectins are also known to play important roles in the immune system. Within the innate immune system, lectins can help mediate the first-line of defense against invading microorganisms. In this study, several kinds of lectins were found to be highly induced in sevenband grouper brain tissue by NNV infection, such as C-type lectins (CTLs), galectins, fucolectin, and mannose-binding lectin. CTLs are the most well studied lectins. They can promote antibacterial and antiviral immune defense [54]. Many CTLs have been identified in teleost but the exact function of CTLs remains unclear.
Hundreds of interferon stimulated genes (ISGs) were transcribed in sevenband grouper brain tissue during NNV infection. Interferon induced protein 44 (IFI44) was expressed the most. IFI44 is an interferon-alpha inducible protein associated with infection of several viruses. Power et al. have demonstrated that IFI44 can inhibit HIV-1 replication in vitro by suppressing HIV-1 LTR promoter activity [55]. Carlton-Smith and Elliott have screened ISGs related to Bunyamwera orthbunyavirus replication using nonstructural (NSs) protein knock out virus. One of these ISGs that have inhibitory activity is found to be IFI44 [56]. Whether protein B2 of NNV has roles in virus replication and its relationship with ISGs such as IFI44 merit further study.
HSPs are one of the most phylogenetically conserved classes of proteins with critical roles in maintaining cellular homeostasis and protecting cells from various stresses [57]. Ironically, HSP70 has roles to suppress some virus infections, and support their replication in other viruses [58]. In this study, HSP30 was the highly induced gene instead of HSP70. HSP30 has also been reported to be the most induced gene in the NNV infected Asian seabass epithelial cell [33]. However, the function of HSP30 in NNV infection remains unclear.
Krasnov et al. previously reported the effects of NNV on gene expression in Atlantic cod brain using a microarray [59]. Compared to our study, a number of genes show a similar up-regulation result in the study, such as Caspase, Cathepsins, IRF, Radical S-adenosyl methionine domain-containing protein, tripartite motif-containing protein (TRIM) and so on. However, a lot of novel genes (e.g., NK-Lysin, Ubiquitin-like protein 4, Granzyme A, etc.) were identified from our RNA-Seq result because a microarray can only evaluate the genes on a chip.
Our findings are preliminary based on the small scale of the study and the results have not yet been confirmed by an independent technique such as quantitative polymerase chain reaction (qPCR). In future studies, it will be necessary to confirm the expression level of genes and to characterize the function of genes that are highly involved in NNV infection.

Conclusions
In conclusion, to the best of our knowledge, this is the first study reporting the transcriptome of brain tissue of NNV-infected sevenband grouper. In this study, we obtained the transcriptome of sevenband grouper. A total of 104,348 transcripts were obtained, including 628 DEGs between NNV infected and non-infected sevenband grouper. A large number of differential expressed genes relevant to immune response were identified as well as several candidate genes (CCL2, Cathepsins, Lectins, Hsp30, and Interferon-induced protein 44) that were intensely induced by NNV. Their functions in sevenband grouper against NNV infection merit further study. The acquired data from such transcriptome analysis provide valuable information for future functional genes related to NNV infection and VNN outbreak.