Gene expression profiles alteration after infection of virus, bacteria, and parasite in the Olive flounder (Paralichthys olivaceus)

Olive flounder (Paralichthys olivaceus) is one of economically valuable fish species in the East Asia. In comparison with its economic importance, available genomic information of the olive flounder is very limited. The mass mortality caused by variety of pathogens (virus, bacteria and parasites) is main problem in aquaculture industry, including in olive flounder culture. In this study, we carried out transcriptome analysis using the olive flounder gill tissues after infection of three types of pathogens (Virus; Viral hemorrhagic septicemia virus, Bacteria; Streptococcus parauberis, and Parasite; Miamiensis avidus), respectively. As a result, we identified total 12,415 differentially expressed genes (DEG) from viral infection, 1,754 from bacterial infection, and 795 from parasite infection, respectively. To investigate the effects of pathogenic infection on immune response, we analyzed Gene ontology (GO) enrichment analysis with DEGs and sorted immune-related GO terms per three pathogen groups. Especially, we verified various GO terms, and genes in these terms showed down-regulated expression pattern. In addition, we identified 67 common genes (10 up-regulated and 57 down-regulated) present in three pathogen infection groups. Our goals are to provide plenty of genomic knowledge about olive flounder transcripts for further research and report genes, which were changed in their expression after specific pathogen infection.

viruses, viral hemorrhagic septicemia virus (VHSV) is affiliated to Novirhabdovirus genus, which is a member of the Rhabdoviridae family 4 . The six gene were contained in the VHSV genome of about 11 K bases and each of them coded nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G), nonstructural viral protein (NV), and RNA polymerase (L) in the following order 3'-N-P-M-G-NV-L-5' 4 . Infection of VHSV results in contagious viral hemorrhagic septicemia (VHS) in diverse fish species regardless of their inhabitation; seawater or freshwater 5 . In East Asia, a lot of infection cases into olive flounder have been reported steadily, since VHSV was detected in middle of 1990s [6][7][8][9] . A variety of scuticociliates have been reported as cause of scuticociliatosis in marine species including turbot, guppy, and southern bluefin tuna [10][11][12] . In olive flounder, disease has been reported to be causing from various scuticociliates; Uronema marinum, Pseudocohnilembus persalinus, Philasterides dicentrarchi, Miamiensis avidus [13][14][15][16] . Interestingly, judging from infection experiments using various scuticociliates plus identification outcome of 8 isolates acquired from olive flounders with symptom of ulcers and haemorrhages, Miamiensis avidus was suggested as the major aetiologic agent of scuticociliatosis because of high pathogenicity and mortality rate compared with other scuticociliates 14,17 .
The main issue of aquaculture industry is to reduce economic loss by preventing mortality of fish from various pathogens. A large number of immunologic studies have been proceeded about various immune-related gens against pathogen infection 3,[24][25][26][27] . A huge quantity of genomic information from next generation sequencing (NGS) technique has been gradually increasing for the last few years, indicating that researchers could approach more comprehensive understanding view about genome of organisms than when they research a single gene level. With development of wide-sized analysis methods, it is not difficult to figure out change of gene expression level after any chemical treatment or environmental change. Recently, studies to identify large-scale genes were conducted in the olive flounder genome for researches about vaccine, gonadal development, and sex determination [28][29][30] . In particular, characterizing of immune-related genes was reported in olive flounder spleen tissue 31 . A lot of studies reported earlier were focused on gene expression analysis of single pathogen and specifically defined the expression pattern of limited genes [32][33][34][35][36] . Further, infection by two or more pathogens were reported in the olive flounder genome 37,38 . In order to solve these problem, we need plentiful genomic information to respond rapidly to multiple infection of pathogens. However, researches, which were comprehensively analysed about change of gene expression pattern by different type of pathogens, have not been reported in the olive flounder genome, so far.
In this research, we identified differentially expressed genes (DEGs) by transcriptome analysis and conducted gene ontology (GO) analysis with genes identified. Then, we tried to find important genes which showed consistently meaningful expression change in the results of three infection experiments. As a result, we determined 10 up-regulated genes and 57 down-regulated genes in common after infection of three pathogens. We aimed to provide essential genome information which is related with pathogen infection and explore the various consequences related to differential infections and find out the common strategies against specific candidates involved in disease progression in natural habitat of aquaculture.

Results
Statistical summary of transcriptome analysis. To profile gene expression after infection of three pathogens (VHSV, Streptococcus parauberis, and Miamiensis avidus), transcriptome analysis was conducted using gill tissues of olive flounders, respectively. We prepared twelve olive flounders (three un-infected individuals as control, three virus-infected, three bacteria-infected, and three parasite-infected individuals) to raise confidence. To gain the sufficient number of transcripts, twelve independent RNA samples acquired from normal and pathogen-infected olive flounder gill tissues were employed for construction of cDNA library. Then, these cDNA libraries were sequenced using Illumina HiSeq2500, generating the numbers of approximately 78.4 million, 65.2 million, and 45.7 million raw reads from three control samples, 44.0 million, 56.9 million, and 62.0 million raw reads from bacteria-infected samples, 41.2 million, 62.4 million, and 40.0 million raw reads from virus-infected samples, 39.6 million, 41.6 million, and 53.1 million raw reads from parasite-infected samples, respectively (Table 1). After trimming of low-quality reads and adaptor sequences, the number of clean reads acquired from control samples were average 62.3 million reads from control samples, average 53.5 million reads from bacteria-infected samples, average 48.0 million reads from virus-infected samples, and average 44.3 million reads from parasite-infected samples, respectively. Then, we checked gene coverage whether the reads that we acquired are sufficient for quantitative gene expression analysis ( Supplementary Fig. 1). The clean reads were assembled into 120,880 transcript sequences acquired from transcriptome analysis, we identified total 40,100 genes involving novel 19,560 genes from transcript sequences using InterProScan database and non-redundant protein database in the NCBI ( Table 2).

Identification of DEGs after viral infection.
To figure out the effects of external pathogen for gene expression, we sorted out genes which showed expressional change after pathogens infection having p-value of <0.05 when compared with control sample. As shown in Table 3 and Fig. 1, the largest numbers of gene expression change were shown in VHSV infection group; total 12,415 DEGs were identified from transcriptome analysis. We showed information of DEGs derived from viral infection in Supplementary  To explore the functional enrichment of these DEGs, we performed GO enrichment analysis using DAVID tool 39 Table 2).

Gene expression change pattern after bacterial infection. After analyzing DEGs and their functional
prediction following viral infection, we identified the DEGs with p-value of <0.05 after infection with bacteria (Streptococcus parauberis) and parasite (Miamiensis avidus) ( Table 3 and Table 2).

Gene expression change pattern after parasite infection. Last, we checked how infection of
Miamiensis avidus affected gene expression pattern in olive flounder genome and selected genes which showed expression change (Supplementary Table 2). We identified 795 DEGs caused by infection of Miamiensis avidus; Description Samples U-I 1 U-I 2 U-I 3 B-I 1 B-I 2 B-I 3 V-I 1 V-I 2 V-I 3 P-I 1 P-I 2 P-I 3

Reads
Number  Distributional pattern of total DEGs acquired from transcriptome analysis. After comparison of gene expression level among twelve transcriptome analysis, we identified DEGs after pathogen infection. Then, we focused on selection of genes showing expression change pattern after three types of pathogen infection in common. As a result, we summarized 10 up-regulated genes and 57 down-regulated genes, respectively (Fig. 2). We analyzed these 67 DEGs to identify their gene symbol correctly using their sequences in non-redundant protein database of the NCBI database, and 37 DEGs were annotated (Table 4). We showed the rest of 30 unannotated DEGs in Supplementary Table 3.

Discussion
With development of sequencing technique, numerous genomic researches have been reported to understand infection results by virus 28,31,40 , bactria 26,41,42 , and parasite 43 in the olive flounder genome. A fundamental way to overcome disease outbreak from external pathogens is to approach from their genome level. It is essential to expand quantity of genomic information in pursuance of biological research about any target. Investigation of overall gene expression change after pathogen infection would provide clues of cause of biological damage. Gene regulation is essential for viruses, prokaryotes and eukaryotes as it increases the versatility and adaptability of an organism by allowing the cell to express protein when needed. Phylogenetic diversity of pathogens (virus, bacteria and parasite) is also responsible for differential expression of genes in diseases. Each individual pathogen causes disease in a different way, which makes it challenging to understand the basic biology of infection. In this study, we understood the relation between three types of pathogen infection and differential gene expression in the olive flounder genome through transcriptome analysis, respectively. The diverse pathogens used in this study, carry specific antigenic variations, which refers to the specific mechanism by which an infectious agent infect the fish and progress the disease. Transcriptome analysis help us to understand the progression of disease in fish through pathogen infection based on diversity of pathogen (virus, bacteria and parasite). This study shows differentially expressed genes were up-and down-regulated   at different extend in fish tissue. Interestingly, virus and bacteria have more down-regulated genes while parasite have more up-regulated genes. This data signifies the fact that fish immune system interacts with bacteria and virus with the same strategy, while with parasite different due to difference in mode of infections between them. For efficient prevention against pathogen, it is important to understand which genes were activated/repressed after pathogen infection because their expression change means variation of metabolic system in body. In this study, we identified total number of 12,415 in VHSV infection group, 1,754 in Streptococcus parauberis infection group, and 795 DEGs from Miamiensis avidus infection group, respectively (Table 3). Given the difference in the number of DEGs among three pathogen groups, these results seemed that virus had the most impact on gene expression mechanism in the olive flounder genome among three pathogens. Interestingly, 11,051 DEGs (89% of all DEGs) showed down-regulated pattern after viral infection. This phenomenon that global gene expression was decreased by viral infection must cause pathogenic disease by affecting immune-related gene expression level, finally leads to death. This view was supported by our findings (Supplementary Table 2), which showed expression decrease pattern of all genes in the immune-related GO terms. Especially, GO terms in viral infection group showed that all genes tend to be down-regulated after pathogen infection, indicating loss of resistance against pathogens by down-regulating the expression of immune-related genes. Like this situation, functional information of genes acquired from GO enrichment could help researcher to figure out critical biological pathway against any external factor.
The immune mechanism in fishes is composed of a set of cellular and humoral system and divided into innate (inherit), and adaptive (acquired) substances. The understanding of fish immune system structure and function is essential for the development of new technologies and products to improve productivity. The transcriptome analysis bring exposure to basic difference in expression profile of all pathogens in the host. These differences were due to nature of parasite, their mode of infection, antigenic variations and many other factors. Additionally, along with all above differences, disease progressed in host due to external surface variation of pathogens (viral, parasite and bacterial) and their appropriate recognition by host immune systems for making the basis to initiate microbial clearance 44,45 . Disease research requires the knowledge of important key factors like method of avoiding host immune surveillance, antigenic variations, subversion of immune responses through phagocyte and inhibition of cytokines and chemokines in common with pathogen infections (viral, bacterial and parasite) 44,45 . On the other hand, the disease progression was different in accordance to type of pathogens. In case of viral infection, understanding of complement inhibition and blockade of cellular immunity is the most important, while parasite and bacterial infections required knowledge and research of innate pathway and acquired immunity 44,45 . Our analysis indicates the complexity and difference of expression profile could be due to all the above reasons. In addition, important basis of fish vaccine is depending on innate and adaptive immunity 46 . There were many vaccine types which depend upon antigens, live microorganisms or specific DNA segment of pathogens or polyvalent vaccines. All above vaccines required complete knowledge of pathogenicity and deep research of efficacy 47 . Our study clearly indicates about various immune and antigenic genes which can be chosen for pathways analysis and use for therapeutic agents or as some vaccine candidates (Supplementary Table 2).
Despite of their different infection pathway, we wondered common DEGs which were affected from infection. As shown in Fig. 2, total 67 DEGs (10 up-regulation and 57 down-regulation) were identified in common in three pathogens, and 1 of 10 up-regulated and 36 of 57 DEGs were annotated by genomic database (Table 4; , 408 genes (bacteria), and 561 genes (parasites). These genes were specific for each pathogen so can be used as candidate genes for vaccination or therapeutic agents. In case of the down-regulation of 10168 gene (virus), 177 gene (bacteria) and 60 gene (parasites) in infection of fish, these genes were specific for specific pathogens so can be used as a diagnosis marker for specific pathogens. C-X9-C Motif Containing 4, COL1A1; Collagen Type I Alpha 1 Chain, COL1A2; Collagen Type I Alpha 2 Chain, SLC14A2; Solute Carrier Family 14 Member 2) which showed highest expression change (down regulation) with fold change (log 2 ) of ' <−2.0' on at least two infection groups. In this study, we investigated all the above candidate genes and found their role in disease progression. We have listed these genes and their role in below headings. ANPEP, called as Gene Aminopeptidase N (APN), is metallopeptidase that exerts strong influence on various immune response mechanisms. For example, APN has been known to cause decomposition of cytokines and peptides used by neurons [48][49][50] , and acts as receptor for viruses 51,52 . In addition, relation between expression level of APN and stimulated T-cell was reported 53 . Recently, it was suggested that APN controlled the balance of innate immune and adaptive immune by regulating TLR4 signal transduction pathway in myeloid cells 54 .
BGLAP, also known as Osteocalcin, is a noncollagenous protein, mainly found in bone, which needs vitamin K for its synthesis. This protein was thought to play a role in calcium ion homeostasis and used as biological marker for bone formation 55 . In addition, it concerns in endocrine regulation, especially in digestive system, by stimulating release of insulin hormone from β-cell of the pancreas and adiponectin hormone from fat cells, respectively 56 . As well as these function, it has been reported to take a role in promotion of energy availability and sexual maturation of male by stimulation of testosterone biosynthesis 57,58 .
CMC4, called as MTCP1, has been mainly reported to be related in various diseases. It was reported that MTCP1 gene affected T-cell homeostasis prior to process of leukemogenesis in transgenic mice 59 . Although the function of this gene has been entirely discovered, regulation error of MTCP1 gene affected on cell survival and cell growth 60 . Besides, this gene was known to be related in the pathogenesis of a subset of T-cell lymphoproliferative diseases 61,62 .
COL1A1 and COL1A2 encodes the pro-α1 chain and the pro-α2 chain protein, respectively. The type I collagen, which is comprised by two pro-α1 chains and one pro-α2 chain, plays a role in reinforcement and support in most of all connective tissues such as bone, cartilage, skin, and tendon, and offers those tissues rigidity and elasticity 63,64 . This protein was reported to stimulate expression of pro-inflammatory cytokines and professional phagocytes in teleost fish gilthead seabream 65,66 . In addition, it has been reported that receptor-mediated interaction which is formed in between cells and collagen molecule might affect in wound healing, inflammatory, and immune response by activating various factors such as cytokines, growth factor, and matrix metalloprotease 63,[66][67][68][69][70] .
SLC14A2, also is known Urea transporter 2 (HUT2), is important gene involved in urea transport and play role in physiology. In mammals, two types of urea transporter (SLC14A1 and SLC14A2) has been reported 71 and were regulated by vasopressin hormone 72 . The kidney uses urea to maintain the appropriate concentration and volume of blood. Without control of these proteins, organism would result in extreme damage in urinary system. Besides, a previous study has reported that genetic variation including nucleotide change is known to significantly influence blood pressure (BP) and metabolism syndrome 73,74 .
As shown in Supplementary Table 2, immune-related DEGs were revealed as results of three pathogens infection. Infection of pathogens caused activation of immune system to respond to invasion of harmful external elements, indicated that change of expression level of immune-related genes. The down-regulation of gene could sequentially influence on expression of various molecules positioned in down-streams in metabolic pathway. In this view, CD13, which is one of down-regulated genes by infection of three pathogens on common, was reported to inactivate interleukin 8 49 . Representative function of this cytokine is to induce migration of neutrophils and granulocytes toward infection site. In addition, absence of CD13 considerably improves cross-presentation of soluble antigen via regulation of receptor-mediated uptake 75 . Thus, decrease of CD13 expression consequentially might activate immune response in the olive flounder. MTCP1 gene induced malignant T-cell transformation 59 and was related in the leukemogenic process of mature T-cell proliferation 61 . This gene was thought to maintain balance of immune response by T-cell. The innate immune system mediates the initial inflammatory response by pathogen infection or injury. For rapid response against external pathogens, infected cells secrete various cytokines to induce effector cells and complements. The type I collagen, which is comprised by proteins coded from COL1A1 and COL1A2, was involved in the expression of pro-inflammatory cytokines in the innate immune system. However, two genes (COL1A1 and COL1A2) showed decreasing expression pattern after infection in our study. Given sampling period (7 days from infection) of olive flounders for this study, it might be explained that the adaptive immune response was activated in the olive flounder genome. It is hard to understand comprehensively about immune system of fish genome. However, our results were expected to contribute for further study by extend of genomic knowledge in the olive flounder.
In conclusion, this study is helpful in understanding infection of the diversified pathogens (antigenic variation) and their role in disease progression in the olive flounder. The differentially expressed genes identified from transcriptome analysis using three types of pathogens could be useful to study the basic diagnosis and therapeutic mechanisms, and offer opportunities for designing the appropriate vaccines or drug targets for pathogen specific candidate genes. Because of lack of genomic information or using one external infection factor, previous studies have been limited to understand global expression pattern of whole genes in the olive flounder genome. We hope that this research would contribute to achieve great outcome in various biological field.

Materials and Methods
Ethical statement. All experiments with the olive flounders in this study were carried out in accordance with the guidelines and regulation approved by Ethical Committee of Pukyong National University.
Preparation of olive flounder gill tissues. Gill tissues from twelve olive flounder (BW = ~50 g, n = 3/ group) including healthy and infected fish with each pathogen were used for this study. Briefly, healthy fish (non-challenged), sampled fish at 7 days post challenge (dpc) with S. parauberis at 5.06 × 10 3 CFU/fish in 1/3 seawater of 21 °C, sampled fish at 7 dpc with VHSV at 106 PFU/fish in 1/3 seawater of 19  Construction of cDNA libraries for transcriptome analysis. Building of transcriptome libraries were conducted by Illumina's TruSeq RNA protocol, and 1-2 μg of total RNA were used in each samples. AMPure XP beads (BECKMAN COULTER) and Ambion Fragmentation Reagents kit (Ambion, Austin, TX) were used for extraction of Poly(A)+ RNA and their fragment, respectively. As the following steps, cDNA synthesis, end-repair, A-base addition, and ligation of the Illumina indexed adapters were carried out according to Illumina's protocol. The size-selected 250-300 bp cDNA fragments were loaded on a 3% Nusieve 3:1 (Lonza) agarose gel for libraries. The cDNA fragments were recovered using QIAEX II gel extraction reagents (Qiagen), and amplified using Phusion DNA polymerase (New England Biolabs) for 14 PCR cycles. The amplified libraries were purified by AMPure XP beads, their concentration and product sizes were assessed on an Agilent 2100 Bioanalyzer. Sequencing of paired-end libraries were conducted with the Illumina HiSeq2500, (2 × 100 nucleotide read length).
Transcriptome analysis and differential gene expression. Transcriptome analysis were carried out with the RNAseq Tuxedo protocol. Mapping of Sequences were conducted against the Olive flounder draft genome (Submitted at present) using TopHat v2.0.9 with default options for paired-end sequences. Transcripts expression were estimated using the Cufflinks program v2.1.1. Total sequencing reads were subjected to preprocessing as follows: adapter trimming was performed using cutadapt with default parameters, and quality trimming (Q30) was performed using FastQC with default parameters. Processed reads were mapped to the Olive flounder draft genome (Submitted at present) using tophat and cufflink with default parameters 76 . The differential analysis was performed using Cuffdiff 76 using default parameters. Further, the FPKM values from Cuffdiff were normalized and quantitated using R Package Tag Count Comparison (TCC) 77 to determine statistical significance (e.g., P values) and differential expression (e.g., fold changes.). Through these statistics analysis, we sorted DEGs having p < 0.05 and showed them as results.
Gene ontology analysis. DEG set for GO analysis was acquired from transcriptome analysis. DEGs were annotated from InterProScan database and non-redundant protein database in the NCBI. DAVID and uniprot tool were used for exploring the functional enrichment of these DEGs and sorting specially out the immune-related GO terms with p-value of <0.05, respectively.