Distinct herpesvirus resistances and immune responses of three gynogenetic clones of gibel carp revealed by comprehensive transcriptomes

Background Gibel carp is an important aquaculture species in China, and a herpesvirus, called as Carassius auratus herpesvirus (CaHV), has hampered the aquaculture development. Diverse gynogenetic clones of gibel carp have been identified or created, and some of them have been used as aquaculture varieties, but their resistances to herpesvirus and the underlying mechanism remain unknown. Results To reveal their susceptibility differences, we firstly performed herpesvirus challenge experiments in three gynogenetic clones of gibel carp, including the leading variety clone A+, candidate variety clone F and wild clone H. Three clones showed distinct resistances to CaHV. Moreover, 8772, 8679 and 10,982 differentially expressed unigenes (DEUs) were identified from comparative transcriptomes between diseased individuals and control individuals of clone A+, F and H, respectively. Comprehensive analysis of the shared DEUs in all three clones displayed common defense pathways to the herpesvirus infection, activating IFN system and suppressing complements. KEGG pathway analysis of specifically changed DEUs in respective clones revealed distinct immune responses to the herpesvirus infection. The DEU numbers identified from clone H in KEGG immune-related pathways, such as “chemokine signaling pathway”, “Toll-like receptor signaling pathway” and others, were remarkably much more than those from clone A+ and F. Several IFN-related genes, including Mx1, viperin, PKR and others, showed higher increases in the resistant clone H than that in the others. IFNphi3, IFI44-like and Gig2 displayed the highest expression in clone F and IRF1 uniquely increased in susceptible clone A+. In contrast to strong immune defense in resistant clone H, susceptible clone A+ showed remarkable up-regulation of genes related to apoptosis or death, indicating that clone A+ failed to resist virus offensive and evidently induced apoptosis or death. Conclusions Our study is the first attempt to screen distinct resistances and immune responses of three gynogenetic gibel carp clones to herpesvirus infection by comprehensive transcriptomes. These differential DEUs, immune-related pathways and IFN system genes identified from susceptible and resistant clones will be beneficial to marker-assisted selection (MAS) breeding or molecular module-based resistance breeding in gibel carp. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3945-6) contains supplementary material, which is available to authorized users.


Background
Gibel carp, also known as silver crucian carp or Prussian carp, has been recognized as a subspecies Carassius auratus gibelio of crucian carp (C. auratus) [1,2], and currently as a separate species C. gibelio owing to its polyploidization and special multiple reproductive modes [3][4][5][6][7][8]. Gibel carp was found to be able to reproduce by unisexual gynogenesis [1,2], bisexual reproduction [4,[9][10][11], hybrid-similar development mode [8] or even androgenesis [12] in response to sperm from different species or gibel carp clones. Since the application of allfemale gibel carp produced by heterologous sperm gynogenesis (termed allogynogenesis) to activate embryo development in the early 1980s [2], gibel carp has become a very important aquaculture species in China and the annual production capacity of crucian carp has increased to 2,912,258 tons in 2015 [13,14]. In aquaculture, gibel carp seeds are generally produced by allogynogenesis that maintains variety purity and high seed survival rate [14]. The homozygous genetic background of each variety has been confirmed by microsatellites, AFLP profiles and transferrin alleles [12]. Recently, the culture industry has suffered enormous economic losses in main culture areas of Jiangsu province due to the epizootic disease caused by crucian carp herpesvirus (CaHV) [15]. CaHV shows high homology to a variant cyprinid herpesvirus 2 (CyHV-2) strain SY-C1 [16] and induces severe mortality. CyHV-2, also named herpesviral haematopoietic necrosis (HVHN), was initially isolated from goldfish (C. auratus) in Japan [17].
Diverse local populations or various gynogenetic clones of gibel carp have been identified in different natural regions of Eurasian continent by karyotypes [18][19][20][21][22][23], serum transferrin phenotypes [24,25], RAPD (random amplification polymorphism of DNA) and SCAR markers [9,20,26], microsatellite markers [27][28][29][30][31][32], transferrin allele polymorphism [25,[33][34][35] or mtDNA control region sequences [25,31,32,[35][36][37][38]. Significantly, some clone-specific molecular markers have been isolated not only for genetic resource identification, but also for marker-assisted selection breeding. In our lab, a series of gibel carp clones, named A, B, C, D, E et al., were successively discriminated from Shuangfeng reservoir, Dongting lake, Pengze lake and other natural regions. Several improved varieties, such as high dorsal allogynogenetic gibel carp (clone D) and allogynogenetic gibel carp "CAS III" (clone A + ), have been successfully bred in the past 30 years [4,10,14]. Currently, the improved variety clone A + is the most popularly cultured varieties in China and the culture scale has accounted for about 70% of gibel carp culture owing to its excellent growth performance [4,10,12,14]. An artificial clone F with subgenomic incorporation was obtained by cold treatment of the clone E eggs inseminated with bluntnose black bream (Megalobrama amblycephala Yin) sperm [19,39,40] and had been proliferated by a dozen successive generations of gynogenesis with Xingguo red common carp (Cyprinus carpio) sperm stimulation, showing rapid growth and disease resistance. Therefore, it is essential to evaluate the susceptibilities of the leading variety clone A + and candidate variety clone F challenged with CaHV. The resistance or tolerance ability of fish to pathogen is determined by genetic factors [41,42], and several disease-resistant varieties were successfully bred in rainbow trout (Oncorhynchus mykiss) [43,44] and Japanese flounder (Paralichthys olivaceus) [45] through selective breeding. Thus, it is an urgent need of gibel carp culture industry to screen resistant gynogenetic clones as core breeding populations to breed novel variety with enhancing CaHV resistance.
Disease resistance is a complex trait that involves various biochemical processes. Host immune responses play central roles in defensing virus attack and involve in innate and adaptive immune systems, such as pathogen recognition receptors, interferon (IFN)-mediated antiviral response, antigen presentation, inflammatory regulators, immune effectors, and so on [46][47][48]. A lot of IFN system genes have been identified from cultured Carassius auratus blastula (CAB) cells induced by UV-inactivated grass carp hemorrhage virus (GCHV) [47,[49][50][51][52][53][54][55][56][57][58][59][60][61][62][63], and their expression regulation and antiviral mechanisms have been revealed in vitro [50,52,54,56,64]. Owing to the complexity of fish antiviral immune responses, it is necessary to gain insights into the nature of antiviral host reactions. RNA-Seq has been proved to be an effective technique to find massive immune-related genes and to better understand the complex interactions between virus and host [65][66][67][68][69]. In this study, we firstly evaluated the susceptibilities of clone A + , clone F and a wild clone H against to CaHV challenge. Based on transcriptome analysis, the common defenses and the distinct immune responses among three clones were revealed. Finally, DEUs upregulated or down-regulated in all three clones or specifically in resistant or susceptible clones were analyzed, and full-scale of expression changes of IFN system genes were investigated. This study provides transcriptomic basis for the mechanism of CaHV resistance and will be beneficial to disease-resistance breeding of gibel carp.

Results
Different resistance of three gibel carp clones in response to herpesvirus infection Three gynogenetic clones of gibel carp, clone A + , clone F and wild clone H were selected to perform herpesvirus challenge experiments by isolated CaHV. Clones A + and H have spindle-shaped body type, while clone F has a bulge on the anterior back just after the head. Clone A + is silver-black in body color, while clone F is silver-white and clone H is yellow-black (Fig. 1a). Moreover, three clones can be discriminated by different transferrin phenotype pattern [25,34,35], and intraclonal homogeneity and interclonal heterogeneity were obvious in these clones (Fig. 1b).
After CaHV infection, the diseased fishes first showed sub-clinical symptoms, including lethargy, anorexia and body color deepened. As the disease progresses, they exhibited bleeding at the base of fins and on abdomen, observed pale gills, internal organ hemorrhaging, pink ascites in abdominal cavity and swollen spleen and kidneys (Fig. 1c). Cumulative mortalities resulted from three independent repetitive experiments were statistically analyzed in these clones. As shown in Fig. 1d, clone A + is most susceptible and clone H is most resistant (p < 0.01). The first death of clone A + occurred at 4 days post injection (dpi) and the overwhelming majority (98.89%) died at 7 dpi. In contrast, the death of clone H started at 6 dpi and about half of individuals (51.11%) died at 14 dpi. Clone F was moderately resistant and 86.67% infected individuals died at 13 dpi. Viral loads were evaluated by real-time PCR. The injected fish from clone H had an average viral load of 10 1.88 and 10 2.16 particles/ng DNA at 3 dpi and 5 dpi respectively, markedly less than those in clone A + and F which had an average viral load of 10 3.29 and 10 2.72 at 3 dpi, and 10 4.49 and 10 3.04 at 5 dpi respectively. Histopathological examination was also performed. In comparison with normal tissue of WT, CaHV infection resulted in severe necrotic lesion, serious vacuolization and hypertrophied nuclei with karyorrhexis in head-kidney (Fig. 1e). These results indicate that clone H possesses stronger resistance to herpesvirus infection.

De novo assembly and functional annotation of headkidney transcriptomes
Head-kidney is the main targeting organ of herpesvirus infection and replication [70]. In order to obtain global gene expression profile responding to CaHV infection, we performed comparative head-kidney transcriptome analysis of control (c) and early diseased (d) individuals with sub-clinical symptoms (an average viral load of 10 6 particles/ng DNA) from clone A + , F and H, respectively (Additional file 1: Figure S1). Three biological replicates were carried out. By using Illumina HiSeq™ 4000 platform, a total of 121.45 Gb data were generated from all transcriptome libraries constructed in this study. After removing the reads with low quality or adaptors, 809,509,754 clean reads were obtained and assembled into a total of 192,369 unigenes with an average size of 1342 bp and a N50 value of 2639 bp ( Fig. 2a; Additional Transferrin phenotype patterns of three clones with three replications. c Symptoms in diseased fish. Hyperemia at the base of the fins and on the abdomen (black arrows), bleeding gills and internal organ hemorrhaging. d Cumulative mortality after CaHV infection. The values are the mean ± SEM from three replicate tanks. Asterisks (*) indicate significant differences (P ≤ 0.05) between clone A + or clone F and clone H. e Histopathological photographs of head-kidney from normal and diseased gibel carp (HE). Healthy fish showed no pathological changes, while CaHV infected fish showed that the head-kidney appeared necrotic lesions (asterisk), serious vacuolization (black triangle), and hypertrophied nuclei with karyorrhexis (blue triangle). Scale bars = 100 um file 2: Table S1, Additional file 3: Table S2) by Trinity [71] and Tgicl [72]. Subsequently, these unigenes were blasted in seven public databases, including National Center for Biotechnology Information (NCBI) nonredundant protein (NR), NCBI non-redundant nucleotide (NT), Gene Ontology (GO), Clusters of Orthologous Groups of proteins (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), InterPro and Swiss-Prot database, and overall 164,017 (85.26%) unigenes were identified (Fig. 2a). Of these, 36,509 unigenes were classified into 62 GO terms (Fig. 2a, b). "Cellular process", "cell part" or "cell" and "binding" were dominant in the category "biological process", "cellular component" and "molecular function", respectively. 1271 unigenes were mapped to the term "immune system process". To identify the pathways involved in physiological function of head-kidney, 93,604 unigenes were mapped to 307 KEGG pathways (Fig. 2a, c). The pathway with most annotated unigenes was "signal transduction" (19,428 unigenes), followed by "immune system" (12,807 unigenes), "cancers: overview" (12,393 unigenes), "infectious diseases: viral" (11,359 unigenes), "infectious diseases: Bacterial" (9652 unigenes), and so on. Therefore, head-kidney triggers violent immune responses when gibel carp is challenged with the herpesvirus CaHV.

Distinct immune responses to the herpesvirus CaHV infection between susceptible and resistant clones
To elucidate the common defense responses and the molecular mechanisms underlying distinct herpesvirus resistance between susceptible and resistant clones, the head-kidney transcriptome profiles of diseased clone A + (d-A + ) versus control clone A + (c-A + ), diseased clone F (d-F) versus control clone F (c-F) and diseased clone H (d-H) versus control clone H (c-H) were compared. A total number of 8772, 8679 and 10,982 DEUs (probability ≥0.8 and relative change ≥2) were identified between the diseased individuals and control individuals from clone A + , F and H, respectively (Additional file 4: Figure  S2; Additional file 5: Table S3).
GO analysis revealed that the DEUs identified from d-A + vs c-A + , d-F vs c-F and d-H vs c-H were categorized into 54, 53 and 54 GO terms respectively, and similar distribution patterns with subtle differences of DEUs were obtained (Additional file 6: Figure S3). Owing to the strikingly more DEUs (>2000) down-regulated in d-H vs c-H (Additional file 4: Figure S2), the numbers of d-H vs c-H DEUs in main GO terms, such as "cellular process", "single-organism process", "metabolic process", "cell", "cell part", "binding" and "catalytic activity", were greater than those identified in d-A + vs c-A + and d-F vs c-F. Of these, 175, 148 and 202 DEUs identified from d-A + vs c-A + , d-F vs c-F and d-H vs c-H were mapped to "immune system process", respectively. Additionally, a few of DEUs were categorized into several GO terms, which were presented only in one clone or two clones, such as "metallochaperone activity" and "hormone secretion" only in d-H vs c-H, "cell killing" in d-F vs c-F and d-H vs c-H, "biological phase" and "rhythmic process" in d-A + vs c-A + and d-H vs c-H, and so on.
KEGG pathway mapping revealed that these DEUs were involved in about 300 similar pathways in three clones. Among the top 30 pathways, 26 pathways enriched in d-H vs c-H were associated with "immune system", "immune diseases", "infectious diseases", "cardiovascular diseases" or immune-related pathways in "signaling molecules and interaction" and "signal transduction", such as "cytokine-cytokine receptor interaction", "Rap1 signaling pathway" and "NF-kappa B signaling pathway" (Fig. 3a; Additional file 7: Table S4). However, 8 pathways among the top 30 pathways enriched in d-A + vs c-A + , including "cell adhesion molecules (CAMs)", "arachidonic acid metabolism", "ECM-receptor interaction", "axon guidance", "osteoclast differentiation" and so on, did not related to immune responses ( Fig. 3b;  Table S4). The differences of enriched pathways between d-H vs c-H and d-A + vs c-A + suggest that differential immune responses occur in susceptible clone and resistant clone.
To reveal the detailed differences of defense responses between susceptible and resistant clones, the numbers of DEUs in KEGG immune-related or disease-related pathways, including "immune system", "immune diseases", "infectious diseases", "cardiovascular diseases", "cancer", "cytokine-cytokine receptor interaction" and "signal transduction involved in immune" were calculated. As shown in Fig. 3d, the numbers of DEUs in "immune diseases", "cardiovascular diseases" and "cytokine-cytokine receptor interaction" were approximately equal among d-A + vs c-A + , d-F vs c-F and d-H vs c-H. However, in term "immune system", "infectious diseases", "cancer" and "signal transduction involved in immune", obviously more unigenes in clone H changed their expression levels than in clone A + and F. Subsequently, the term "immune system" were further sub-classified into 16 pathways, such as "hematopoietic cell lineage", "complement and coagulation cascades", "antigen processing and presentation", "Toll-like receptor signaling pathway" and so on, and "signal transduction involved in immune" included "Jak-STAT signaling pathway", "Rap1 signaling pathway", "NF-kappa B signaling pathway", "PI3K-Akt signaling pathway" and "TNF signaling pathway" (Fig. 3e). The numbers of DEUs were approximately equal between

d-A + vs c-A + , d-F vs c-F and d-H vs c-H are abbreviations of clone A + diseased fishes (d-A + ) versus clone A + control fishes (c-A + ), and so on
two susceptible clones A + and F in all 21 pathways. In contrast, significantly more DEUs were identified from d-H vs c-H than those in d-A + vs c-A + and d-F vs c-F, except "complement and coagulation cascades", "antigen processing and presentation", "cytosolic DNA-sensing pathway" and "PI3K-Akt signaling pathway". The term "infectious diseases" was also subdivided into 10 pathways in "infectious diseases: bacterial", 6 pathways in "infectious diseases: parasitic" and 7 pathways in "infectious diseases: viral". Compared to the faint differences of DEUs numbers among three clones in "infectious diseases: bacterial" and "infectious diseases: parasitic", the number of DEUs in "infectious diseases: viral" identified from d-H vs c-H were remarkably greater than those in d-A + vs c-A + and d-F vs c-F (Fig. 3f), consistent with the herpesvirus nature of CaHV. All the results suggest that distinct immune responses occur between resistant and susceptible gibel carp clones after the herpesvirus CaHV infection.
Hierarchical clustering of the common DEUs was classified into two distinct clusters (A and B) according to the different changes among three clones (Fig. 4c). Cluster A included the commonly up-regulated DEUs and was sub-divided into five distinct clusters. The DEUs in cluster 1 were highly up-regulated in d-A + vs c-A + and d-F vs c-F, in which the vast majority of the DEUs were annotated as intestinal mucin-2. The DEUs in cluster 2 and 4 showed relatively more increases in resistant clone H than in susceptible clone A + , including IFNphi2, IFI56, IFI58, Mx1, Mx3, viperin, Gig1, IRGE4, IL-6a and E3 ubiquitin-protein ligase TRIM39, which might be important candidate resistant genes to CaHV for further studies. The DEUs in cluster 3 displayed relatively higher increase levels in d-F vs c-F than in d-A + vs c-A + and d-H vs c-H, and the DEUs in cluster 5 showed similar increase levels among three clones after CaHV infection. The commonly down-regulated DEUs were categorized into seven distinct clusters in cluster B, in which cluster 7 included the vast majority of DEUs showing similar a little reduction among three clones after CaHV infection. Cluster 4 contained the DEUs highly down-regulated, such as IRF5, IL6R, CFD, c-lectin and so on. The DEUs in cluster 2, 5 and 6 displayed the most expression decreases in clone A + , F and H respectively. The representative genes in cluster 6 included CCR7, CD59, leukotriene B4 receptor Ltb4r1, and so on. The DEUs in cluster 1 or 3 exhibited more decrease in clone A + and F than clone H or in clone A + and H than clone F, respectively.
Specifically changed DEUs and differential immune pathways of three different clones in response to the herpesvirus CaHV infection  Table S6). The KEGG mapping analysis showed the differential enriched pathways among these DEUs ( Fig. 5; Additional file 8: Table  S5). The DEUs specifically upregulated in susceptible clone A + were chiefly mapped to "intestinal immune network for IgA production", "hematopoietic cell lineage", "platelet activation", "leukocyte transendothelial migration", "PI3K-Akt signaling pathway", "HIF-1 signaling pathway", "cytokine-cytokine receptor interaction" and "Toll-like receptor signaling pathway", while the DEUs specifically upregulated in resistant clone H were mainly mapped to "cytosolic DNA-sensing pathway", "hematopoietic cell lineage", "intestinal immune network for IgA production", "complement and coagulation cascades", "RIG-I-like receptor signaling pathway" and "NF-kappa B signaling pathway". The DEUs specifically upregulated in clone F were mapped to "Toll-like receptor signaling pathway", "NOD-like receptor signaling pathway" and "cytosolic DNA-sensing pathway". Interestingly, there exist remarkable differences in terms of down-regulated DEUs of three gynogenetic clones. Among the top 30 KEGG pathways enriched by down-regulated DEUs specific to resistant clone H, 9 pathways were associated with "immune system". However, only 4 and 5 KEGG pathways enriched by down-regulated DEUs specific to clone F and A + were clustered into "immune system". Besides the shared 3 pathways (e.g., "leukocyte transendothelial migration", "hematopoietic cell lineage" and "natural killer cell mediated cytotoxicity") in three gynogenetic clones and "intestinal immune network for IgA production" shared in clone F and A + , other 2 pathways, "B cell receptor signaling pathway" and "complement and coagulation cascades", were down-regulated in susceptible clone A + , while 5 different pathways (e.g., "T cell receptor signaling pathway", "chemokine signaling pathway", "Fc epsilon RI signaling pathway", "platelet activation" and "Fc gamma R-mediated phagocytosis") were observed to be enriched in resistant clone H. Additionally, 229 and 102 DEUs specifically down-regulated in clone H mapped to "Rap1 signaling pathway" and "cytokine-cytokine receptor interaction", respectively. The differences of enriched pathways reveal a marked differences in immune response of gibel carp three gynogenetic clones.
The numbers of these DEUs in KEGG immune-related or disease-related pathways were also calculated. As shown in Fig. 6a, obviously more DEUs in "immune system", "infectious diseases", "cancer" and "signal transduction involved in immune" were detected in clone H than in clone A + and F. Consistent with the results of all DEUs identified from d-A + vs c-A + , d-F vs c-F and d-H vs c-H (Fig. 3d), the number of DEUs unique in resistant clone was greater than those in susceptible clone A + and F, except "complement and coagulation cascades", "antigen processing and presentation", "cytosolic DNA-sensing pathway" and "PI3K-Akt signaling pathway" (Fig. 6b).

Significantly changed IFN system genes of three different clones in response to the herpesvirus CaHV infection
To globally investigate the expression changes of IFN system genes, we searched them from the de novo transcriptome assembly data (Additional file 3: Table S2). The DEUs annotated as IFN system genes and exhibiting representative expression pattern were selected to perform hierarchical clustering analysis. As shown in Fig. 7, IFN system genes were categorized into three distinct clusters (A, B and C) according to their expression change status. The genes in cluster A and cluster C displayed contrasting profiles of gene expression after CaHV infection. The cluster A was sub-divided into 6 clusters. The DEUs in cluster 1, 4 and 5, including IFNγ-1, IFNγ-2, IFNphi2, IFNphi3, Mx1-like-1, IFI44-like-1, IFI44-like-2, IFI44-like-3 and Gig2, remarkably or relatively increased their expression levels among three clones. Interestingly, IFNphi3, IFI44-like-1, IFI44-like-2, IFI44-like-3 and Gig2 showed the highest increase in clone F after CaHV infection. Cluster 2 and 3 consisted of the resistant-related genes, such as viperin, PKR-1, PKZ, IRF7, Mx1-1, IFI56, Mx-3, Gig1-1, Gig1-2, RIG-I, and so on, which displayed relatively more increases in resistant clone H than in susceptible clone A + or specifically increased only in resistant clone H. IRF1 was very special and uniquely raised its expression only in susceptible clone A + . The cluster C was sub-divided into 4 clusters. The DEUs in cluster 2 and 4 reduced their expression among three clones, consisting of IFNphi4, IFI30, TBK1-1, IRF5-1, IRF5-2 and so on. Cluster 2 included the DEUs displaying relatively more decreases in resistant clone H than in susceptible clone A + , while the genes in cluster 3 showed remarkably more reductions in resistant clone H, such as JAk2-1, TLR2, TLR5, IFITM1, and so on. The genes in cluster B, including IFN1, IRF2, IRF3, IRF6, STAT1, STAT2, interferon alpha/beta receptor IFNAR1, interferon-inducible double stranded RNA dependent inhibitor PRKRIRA, interferon-inducible double stranded RNA dependent activator PRKRA, interferon stimulated exonuclease gene ISG20L2, interferon-related developmental regulator IFRD1, TANK-binding kinase TBK1, melanoma differentiation-associated gene MDA5, myeloid differentiation primary-response protein MyD88, and so on, showed little change after CaHV infection among three clones.

Discussion
Currently, the production capacity of gibel carp contributes to approximate 10% of Chinese freshwater aquaculture production [14]. Owing to the highly intensifying and seriously crowding of monoculture, gibel carp has been suffered by a series of pathogens [73]. Since 2009, an epizootic with acute gill hemorrhages and high mortality has outbroken in cultured gibel carp [74]. The complete genome sequencing of the virus isolated from the tissues of diseased gibel carp showed the CaHV was most closely related to CyHV-2 [15]. CyHV-2 was first reported in 1992 as a pathogen of goldfish in Japan [17,75], then was identified from goldfish in many countries, such as USA [76,77], Australia [78], UK [79] and New Zealand [80]. By using PCR assay, CyHV-2 was first detected in gibel carp in Hungary [81]. Recently, it has been detected in many main culture regions of gibel carp [74] and caused huge economic loss. The experimental infection of indigenous Cyprininae species in Japan, such as ginbuna C. auratus langsdorfii, nagabuna C. auratus buergeri, nigorobuna C. auratus grandoculis and common carp, suggested that Japanese Carassius fish species possess different ability to depress the replication of CyHV-2 [82]. To obtain gynogenetic clones with strong resistance to CaHV, we evaluated the susceptibilities of a leading variety A + in China, a candidate variety F and wild clones in responsive to CaHV. The cumulative mortalities (Fig. 1d) and viral loads in infected fishes both indicate wild clone H has higher resistance to CaHV than clone A + and F and can be used as core breeding populations to breed novel variety with enhancing CaHV resistance.
Consistence with KEGG pathway mapping analysis, a lot of IFN system genes, including IFNγ, IFNphi2, IFI35, IFI44-like, IFI56, IFI58, Mx1, Mx3, MxE, viperin, Gig1, Gig2 and GBP1, were up-regulated in all three clones (Additional file 9: Table S6), indicating that IFNmediated innate immune response is one of major common immune defenses of gibel carp to CaHV. Significantly, most of these IFN system genes, such as viperin, PKR-1, PKZ, IRF7, Mx1-1, IFI56, Mx-3, Gig1-1, Gig1-2, RIG-I, and so on, showed relatively more increases in resistant clone H than in susceptible clone A + or specifically increased their expression only in resistant clone H (Fig. 7). The similar dynamical expression changes of IFN system genes were also observed between resistant and susceptible Atlantic salmon (Salmo salar) families challenged with infectious pancreatic necrosis virus (IPNV) [90,91]. At 1, 5 and 21 dpi, the expression levels of IFN, Mx1 and PKR in headkidney of resistant family were higher than those in susceptible family [90]. In another resistant and susceptible full-sibling families of Atlantic salmon, IFNα and IFNγ increased rapidly at 1 dpi, then dropped to the basal values at 5 dpi in susceptible family, while their expression slightly raised and then maintained in resistant family at 5 dpi which were higher than those in susceptible family [91]. Conversely, a completely different immune reaction was observed in other resistant and susceptible families of Atlantic salmon. The vast majority of innate immune response genes, including IFNα, IFNγ, Mx, ISG15, viperin and Gig2, had higher expression in whole fry of susceptible family than in resistant family [92]. These different dynamical changes of IFN system genes among different Atlantic salmon families might be due to different tissues used (head-kidney in the two formers and whole fry in the latter) or different families [92], which indicated that fish immune response against virus is complex and possesses a species-specific or tissuespecific manner. Further studies are necessary to verify the changes of these important DEUs identified in this study at different stages after infection and to reveal their regulative mechanisms behind differential expression between gibel carp resistant and susceptible clones. And the association between single nucleotide polymorphisms (SNPs) of immune-related genes (IL-10a and MHC class IIB) and resistance to the cyprinid herpesvirus-3 (CyHV-3) had been revealed in common carp [93,94]. Thus, the differential DEUs, especial IFN system genes showing upregulation expression in resistant clone H after CaHV infection, might be a key factor for its stronger resistance to CaHV and their allelic variation related to herpesvirus resistance will be able to use as molecular module markers for disease-resistance breeding in gibel carp.
Inflammatory response is crucial for protection against pathogens, and causes significant tissue damage, which involved in cytokines. Besides IFN, cytokines also include chemokines, interleukin, tumor necrosis factor super-families (TNF), colony stimulating factors and so on [95,96]. IL-6 superfamily is produced in the early stages of infection and involves in diverse immune and neuroendocrine processes, including the regulation of lymphocyte and monocyte differentiation, migration of leukocytes towards the sites of inflammation and chemokine secretion [97]. In this study, the pathway with most annotated unigenes of commonly up-regulated DEUs in KEGG was "cytokine-cytokine receptor interaction". Two members of IL6 superfamily, IL6 and IL11, raised their expression in all three clones (Additional file 9: Table S6). Significantly, IL6 increased more in resistant clone H than in susceptible clone A + . The up-regulation of IL6 was also observed in Japanese pufferfish (Fugu rubripes) [98], rainbow trout [99], Japanese flounder [100], gilthead seabream (Sparus aurata) [101], and European sea bass (Dicentrarchus labrax) [102]. After LPS, poly I:C or pathogen infection, Japanese flounder, rainbow trout and gibel carp IL11 also increased expression after challenge with viral [74,103,104]. Rainbow trout recombinant IL6 protein can promote macrophage growth in culture and induce up-regulation of antimicrobial peptide [105]. The immunoglobulin production is mainly regulated by IL6 [106]. Fugu (Takifugu rubripes) or orange-spotted grouper recombinant IL6 protein can induce the production of IgM [107,108]. In accord with the higher expression of IL6 in resistant clone than in susceptible clone, a lot of DEUs annotated as immunoglobulin were identified from the file including up-regulated DEUs specific to d-H vs c-H (Additional file 9: Table S6). In contrast, IL10 and IL22 were found in the up-regulated DEUs specific to d-A + vs c-A + (Additional file 9: Table S6). IL10 was initially identified as a cytokine synthesis inhibitory factor [109] and had been reported to downregulate IL6 expression [110]. Indian major carp (Labeo rohita) recombinant IL10 protein induced down-regulation of most pro-inflammatory cytokines and up-regulation of natural killer enhancing factors [111]. The higher up-regulation of IL10 in clone A + suggests a weaker immune response. Additionally, IRF1 specifically raised its expression level only in susceptible clone A + after CaHV infection. IRF-1 promotes apoptosis following DNA damage [112]. Consistent with upregulation of IL10 and IRF-1 in susceptible clone A + , a lot of apoptosis or death related genes, such as p53, bcl-2-like, Mcl1b, CD244, dapk3 and FAS, showed more increases or specifically raised their expression in clone A + after exposing to CaHV (Additional file 9: Table S6).
After exposure to CaHV, a lot of DEUs mapped to "complement and coagulation cascades" were downregulated in all three gibel carp clones ( Fig. 4b; Additional file 8: Table S5). The relationship between virus and complements is very complicated [113][114][115].
The complement system plays multiple roles in defending virus evasion, including alerting host to the presence of virus, eliminating invading virus, promoting inflammatory response, clearing apoptotic cell and necrotic cell debris, and modulating innate and adaptive immune response [113][114][115][116][117][118][119][120]. In the meantime, virus has evolved multiple strategies to escape the attack of complements, such as expressing mimics of host complement regulators [113][114][115]. Although the mechanisms of complement activation and evasion have been revealed clearly in mammals [113], the elaborate interactions between complements and viruses in fish are still unknown and the dynamic expression changes of fish components after virus infection are inconsistent. Generally, fish components are up-regulated after virus infection. However, several exceptions indicate complement system might act various immune effecter functions depending on the nature of pathogens or the expressing tissue [116,121]. Complements (e.g., C7 and Df) were downregulated in rainbow trout after viral hemorrhagic septicaemia virus (VHSV) [122]. In zebrafish, the infected individuals showed down-regulation of C3, C8a, C8g, Crpp and Hf in internal organs after VHSV bath challenge, while up-regulation of many complements at the fins [123]. In addition, the expression levels of factor B and C3 were remarkably lower in the resistant group than in susceptible group of Japanese flounder (Paralichthys olivaceus) vaccinated with Streptococcus iniae [124]. In this study, the expression of C3, C4-1, C5-2, C6, C7, D and other complements were down-regulated in all three clones of gibel carp (Additional file 9: Table S6). Moreover, complements and TLR signaling pathways can influence each other in regulating inflammatory responses [125,126]. Recently, Lebel et al. [127] observed increased immune cell activation and higher production of IFN-α in C3-depleted mice treated by papaya mosaic virus (PapMV)-like nanoparticle, consistent with the activation of IFN system and suppression of complement system in gibel carp after CaHV infection.
A large number of DEUs were mapped to pathway "leukocyte transendothelial migration", "T cell receptor signaling pathway" and "B cell receptor signaling pathway". Leukocyte transendothelial migration is a critical step in immune activation. Following immune activation, several chemokines were up-regulated, while chemokine receptors and interleukin receptors were down-regulated their expression after CaHV infection. Chemokines belong to a family of structurally related chemotactic cytokines and regulate migration of monocytes, neutrophils and other effector cells to the sites of tissue infection [128,129]. The similar dynamical expression changes were observed in other fishes. After bacterial infection with Edwardsiella ictaluri, 9 channel catfish (Ictalurus punctatus) chemokines were up-regulated [130]. In large yellow croaker (Pseudosciaena crocea), CXCL12 was significantly up-regulated in many tissues after stimulation [131]. In addition, a lot of CD antigens expressed on leukocytes and other cells relevant to the immune system were found to change their expression after CaHV expression. More than 400 CD molecules had been identified in human and classified into about 50 superfamilies, including immunoglobulin superfamily (IgSF), G-protein coupled receptor superfamily, C-type lectin family, cytokine receptor family, TNF superfamily, TNF receptor superfamily, integrin family, tetraspanin family, Toll-like receptor family, cadherin family, and so on [132,133]. In this case, CD molecules belonging to IgSF (e.g., CD2, CD3, CD4, CD8, CD22, CD79b, CD80, CD86 and CD276), cytokine receptor family (e.g., CCR4, CCR5 and CCR7), TNF/TNFR superfamily (e.g., CD40), integrin family (e.g., CD11 and CD18) and tetraspanin family (e.g., CD9 and CD81) were down-regulated at least in one clone of gibel carp after CaHV infection.
Although the biological functions of CD antigens have been well investigated in mammals, their characterizations in fish remain unclear. Zebrafish (Danio rerio) CD44 and CD154 were significantly up-regulated after stimulating with KLH [134], while CD36 was downregulated during infection of Mycobacterium marinum [135]. In rainbow trout, CD antigens in responsive to pathogens are complicated. After VHSV infection or stimulation of Vibrio bacterin, CD28 expression was not significantly different in control and infected fish, while CD152 was up-regulated in spleenocytes [136], CD80/86 expression raised in leukocytes [137] and CD3, CD8 and CD4 increased expression in liver [138] respectively. Additionally, CD9 and CD63 were down-regulated in gill or head-kidney by VHSV bath challenge, while they both significantly up-regulated in peritoneal cells when the virus was intraperitoneally injected [139]. In rock bream (Oplegnathus fasciatus), the expression dynamical changes of CD200 depended on pathogens or tissues [140]. Therefore, more studies are needed to confirm the changes of complements, chemokines, and CDs in gibel carp clones.

Conclusions
This work is the first report to obtain resistant clone against the herpesvirus CaHV and the first transcriptomic comparison between susceptible and resistant clones in gibel carp. Taken together, gibel carp activates IFN system and suppresses complements and CD antigens to defend CaHV invasion. Resistant clone H triggers stronger immune responses with higher expression of significant IFN system genes, IL6 and immunoglobulins, while susceptible clone A + fails to protect from CaHV infection with more expression of apoptosis or death related genes. The further studies of resistance-relevant genes identified in this study could provide useful information for disease control with effective immune protection and for resistance breeding in gibel carp.

Fish
Six month old gibel carps (Carassius gibelio) were obtained from the GuanQiao Experimental Station, Institute of Hydrobiology, Chinese academy of sciences, which is located in Wuhan, China. The average weights of clone A + , F and H were 85.89 ± 2.13 g, 67.68 ± 2.16 g, 46.96 ± 2.47 g respectively. Apparently healthy individuals were selected to gradually acclimatize in 150-l tanks with aerated water at 24(±1)°C for 2 weeks before infection and fed with commercial feed twice a day. Before infection, 5 individuals of each clone were randomly selected to perform PCR analysis as previously described [74] to confirm CaHV free.

Transferrin phenotype
Sera was collected from blood of 3 individuals of each clone by centrifugation and treated with rivanol to isolate transferrin, which was applied to 10% polyacrylamide gel electrophoresis (PAGE) following the procedure by Li and Gui [35].

CaHV infection and sample collection
The CaHV was isolated from the tissues of naturally diseased gibel carp with acute gill hemorrhages. The CaHV was amplified by injection with tissue filtrate from naturally infected fish into healthy gibel carp, identified by PCR assay and electron microscopic observation, and titrated by real-time PCR analysis as previously described [70,74]. Individuals of each clone were randomly divided into six tanks with 30 fishes per tank. Five tanks of each clone were performed infection with 500 μL CaHV viral suspension (2.915 × 10 8 virus particles) per fish by intraperitoneal injection. The remaining fishes were injected in equal phosphate-buffered saline solution (PBS) as the control group. Three of the infected tanks were used to record mortality and the rest were used to collect samples. After infection, fish were put back and grown at 24(±1)°C. Water was filtered consecutively and changed daily to keep it visibly clean. Dead fish were taken out timely. Head-kidney tissues were collected from control fishes, infected fishes at 1, 3, 5 and 7dpi, and fishes with sub-clinical and clinical symptoms respectively. The experiment was terminated after 28 days post infection when mortality stabilized. All samples were preserved in RNAlater (QIAGEN) and stored at −20°C for nucleic acid extraction.

Histopathology
Head-kidney tissues from 3 individuals of control and diseased fishes with clinical symptoms were fixed in 4% paraformaldehyde over night at 4°C. After dehydrated and embedded, samples were cut into 4 μm sections and stained with Hematoxylin-Eosin as previously described [70].

DNA extraction and quantification
Five infected individuals of each clone at 1, 3, 5 and 7 dpi, and diseased fishes with sub-clinical symptoms were selected to evaluate viral load. Total DNA was extracted using DNA extraction kit (Promega, USA) according to the manufacturer's protocol. Quantification of viral copy numbers were calculated by real-time PCR analysis as previously described [74]. Briefly, a 637-bp helicase gene fragment of CaHV was amplified to serve as the standard for virus quantification. The amplified fragment was purified using a Gel Extraction Kit (OMEGA) and inserted into the pMD18-T plasmid to produce pMD-CaHV (3329 bp). A 10-fold dilution series of pMD-CaHV was used as the standard template of CaHV in the quantitative real-time PCR. The primers and procedure of real-time PCR were designed as described [74].

RNA extraction and RNA-Seq
In order to eliminate the differences in virus susceptibility among individuals, head-kidney tissues from 3 individuals with sub-clinical symptoms and similar average viral load (10 6 particles/ng DNA) of each clone were collected and performed the transcriptome analysis. Total RNAs were isolated using SV Total RNA isolation System (Promega, USA) according to the manufacturer's protocols. The quantity and quality of total RNAs were assessed by Agilent 2100 Bioanalyzer using Agilent RNA 6000 Nano Kit and agarose gel electrophoresis. The subsequent library construction, sequencing and bioinformatics analysis were accomplished by Beijing Genomics Institute (BGI), China. Briefly, total RNAs were digested by Dnase I (NEB) and purified by oligo-dT beads (Invitrogen), then fragmented with Fragment buffer (Ambion). The first strand cDNAs and second strand cDNAs were synthesized successively by First Strand Master Mix and Second Strand Master Mix (Invitrogen). The cDNAs were purified and combined with End Repair Mix. After purified, the end-repaired DNAs were mixed with A-Tailing Mix and combined with Adenylate 3′Ends DNA, Adapter and Ligation Mix. A narrow 300 bp-350 bp size-range DNA fragments were selected and enriched. The final libraries were quantified by realtime qPCR (TaqMan Probe) and the average molecule length determined by using the Agilent 2100 bioanalyzer instrument (Agilent DNA 1000 Reagents). The qualified libraries were amplified on cBot to generate the cluster on the flowcell. Then the amplified flowcells were sequenced pair end on the HiSeq 4000™ System and 150 bp single-end reads were generated.

Sequence assembly and annotation
After removing the reads with low quality or adaptors, de novo assembly was performed with clean reads to produce the unigenes for downstream bioinformatics analysis. Quality control of clean data was performed through drawing base composition chart and quality distribution chart. Then the unigenes were aligned to NT, NR, COG, KEGG and Swiss-Prot database using Blast, GO database using Blast2GO [141], and InterPro database using InterProScan5 [142]. With the GO and KEGG annotation, unigenes were classified according to official classification and the functional enrichment were performed using phyper, a function of R. False discovery rate (FDR) was used to determine the threshold of p value and GO or KEGG terms (FDR ≤ 0.01) were considered significantly enriched.

Differential expression analysis
To quantify the unigenes expression level, clean reads were mapped to unigenes using Bowtie2 [143] and the expression levels were calculated by RSEM (RNASeq by Expectation Maximization) [144]. Then DEUs were detected with NOIseq based on noisy distribution model [145] and the threshold to judge the significant expression difference was performed by "probability ≥ 0.8 and relative change ≥ 2". To identify the pathways that DEUs participate in, KEGG enrichment analysis was performed by phyper. Heatmap was generated using MeV and venn diagram was created using jvenn [146].