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

Abstract

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.

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 all-female 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 up-regulated 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).

Fig. 1

Heterogeneity and different susceptibility to CaHV among three gibel carp clones. a Morphological difference among clone A⁺, F and H. b 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

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 head-kidney 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⁶ 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 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) non-redundant 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.

Fig. 2

Annotation summary and unigene classification of the gibel carp transcriptome. a Overview of unigene annotation. b Histogram of GO classification. 36,509 unigenes were assigned to 62 terms grouped into three main categories: biological process (red bars), cellular component (blue bars) and molecular function (green bars). c Histogram of KEGG classification. 93,604 unigenes were mapped to 42 terms grouped into six main categories: Cellular Processes (red bars), Environmental Information Processing (orange bars), Genetic Information Processing (green bars), Human Diseases (dark green bars), Metabolism (blue bars) and Organismal Systems (purple bars)

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; Additional file 7: 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.

Fig. 3

KEGG pathway enrichment analysis in each clone. Top 30 enriched KEGG pathways of DEUs from d-H vs c-H a, d-A⁺ vs c-A⁺ b and d-F vs c-F c. The x-axis indicates the rich factor of each pathway and y-axis shows pathway. The color and size of dot indicates Q value and the numbers of DEUs assigned to the corresponding pathway respectively. d The numbers of DEUs in immune-related or disease-related KEGG pathways at level 2. e The numbers of DEUs in immune-related KEGG pathways at level 3. f The numbers of DEUs in infectious diseases-related KEGG pathways at level 3. 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

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 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.

Common defense responses of three different clones in response to the herpesvirus CaHV infection

The shared DEUs among d-A⁺ vs c-A⁺, d-F vs c-F and d-H vs c-H were integrated. 5465 up-regulated and 15,843 down-regulated DEUs were obtained. Among these DEUs, 407 and 1596 DEUs were up-regulated or down-regulated in all three gibel carp clones after CaHV infection (Fig. 4a). The pathway with most annotated unigenes of commonly up-regulated DEUs in KEGG was “cytokine-cytokine receptor interaction” which indicates cytokines play a vital role in resistance to CaHV invasion. In addition, “Toll-like receptor signaling pathway”, “cytosolic DNA-sensing pathway” and “Jak-STAT signaling pathway”, were also observed in the top 30 enriched KEGG pathways. In contrast, “complement and coagulation cascades”, “intestinal immune network for IgA production”, “hematopoietic cell lineage”, “leukocyte transendothelial migration” and “chemokine signaling pathway” were found in the top 30 enriched KEGG pathways of commonly down-regulated DEUs (Fig. 4b; Additional file 8: Table S5). The results suggest that CaHV infection result in the activation of IFN system genes and the down-regulation of immune related genes involved in complement system.

Fig. 4

Common DEUs among gibel carp clones. a Venn diagram of up- or down-regulated DEUs among three comparisons: d-A⁺ vs c-A⁺ (green circle), d-F vs c-F (blue circle) and d-H vs c-H (red circle). b Top 30 enriched KEGG pathways of commonly up- or down-regulated DEUs in clone A⁺, F and H. The x-axis indicates the rich factor of each pathway showed at y-axis. The color and size of dot indicates Q value and the number of DEUs assigned to the corresponding pathway respectively. c Heatmap of commonly changed DEUs in three clones with hierarchical clustering calculated by log₂ fold change values and the main clusters were lined out on the right. Up-regulated DEUs are shown in red and down-regulated DEUs in green

Approximately 38.3% (156/407) up-regulated DEUs or 34.5% (551/1596) down-regulated DEUs in all three clones were mapped into KEGG immune-related pathways or disease-related pathways (Additional file 8: Table S5). A lot of IFN system genes, such as IFN (IFNγ, and IFNphi2), interferon-induced protein (IFI35, IFI56 and IFI58), IFN antiviral effectors (Mx1, Mx3 and MxE), IFN-stimulated genes (ISGs) (viperin and Gig1), interferon-induced guanylate-binding protein GBP1, and immunity-related GTPase family IRGE4, interleukins (IL-6a and IL-11), chemokines (CXCL-C1c, CCL-C5a and CXCl1) and tumor necrosis factor superfamily members (tnfrsf1a, tnfrsf5, tnfrsf11b, tnfrsf12a, tnfrsf13b and tnfrsf21) were up-regulated in all three clones. Among the 551 down-regulated DEUs, a lot of unigenes belonged to “complement and coagulation cascades”, such as complement factor properdin (CFP), CFB, CFD, CFI, CFHR3, CFHR-1, C3, C4–1, C5–2, C6, C7, complement C3b/C4b receptor 1 (CR1), or were clustered into “antigen processing and presentation”, such as different alleles of major histocompatibility complex (MHC) class II alpha or beta antigen, MHC class II transactivator (CIITA) isoform X2 and isoform X3. Additionally, the expression of several cluster of differentiation (CD) antigens (CD3, CD4, CD9, CD44, CD59, CD79b and CD209), C-C chemokine receptors (CCR4, CCR5 and CCR7), interleukin receptors (IL1R, IL2Rβ, IL6Rα, IL13Rα2), IFN system genes, including IFNphi4, IFI30, dsRNA signal sensing factor (TLR3 and TLR25) and interferon regulator factor (IRF2, IRF4, IRF5 and IRF8), were decreased after 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

As shown in Fig. 4a, 4681 (1510 up-regulated and 3171 down-regulated), 4676 (1466 up-regulated and 3210 down-regulated) and 6829 (1302 up-regulated and 5527 down-regulated) changed their expression levels in clone A⁺, F and H respectively after the herpesvirus CaHV infection (Fig. 4a; Additional file 9: 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.

Fig. 5

KEGG pathway enrichment analysis of unique DEUs in clone A⁺, H or F, respectively. Top 30 enriched KEGG pathways of specifically up- or down-regulated DEUs in clone A⁺, H or F. The x-axis indicates the rich factor of each pathway. The color and size of dot indicates Q value and the number of DEUs assigned to the corresponding pathway respectively

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).

Fig. 6

DEUs in corresponding KEGG pathways. a Histogram of the numbers of specifically changed DEUs in immune-related or disease-related KEGG pathways at level 2. b Histogram of the numbers of specifically changed DEUs in immune-related KEGG pathways at level 3

Among these DEUs specific to d-H vs c-H, approximately 15.1% (196/1302) up-regulated DEUs and 21.5% (1191/5527) down-regulated DEUs were mapped into KEGG immune-related pathways or disease-related pathways (Additional file 8: Table S5). A serial of IFN system genes, such as interferon-inducible and dsRNA-dependent protein kinase PKR, PKZ, IFI27-like, IRF7, stimulator of interferon MAVS, and interferon-induced protein with tetratricopeptide repeat IFIT, specifically increased their expression after CaHV infection, while several other IFN system genes or immune-related genes, such as Toll-like receptor 5 (TLR5), Janus kinase 2 (JAK2), and interferon induced transmembrane protein 1 (IFITM1), uniquely decreased their expression. In addition, after exposing to CaHV, immunoglobulin (Ig) light chain, Igh, Ighm, Ighv, Igll1, polymeric immunoglobulin receptor-like (pigrl) and a novel immune type receptor with Ig domains specifically raised their expression, while CD antigens (CD8, CD11, CD18, CD22, CD80, CD81, CD247 and CD276), C-C chemokine receptors (CCR3.1, CCR9 and CCR10) and interleukin receptors (IL4R, IL7Rα and IRAK4) uniquely reduced their expression.

Among these DEUs specific to d-F vs c-F, approximately 20.1% (293/1466) up-regulated DEUs and 15.7% (505/3210) down-regulated DEUs were clustered into KEGG immune-related pathways or disease-related pathways (Additional file 8: Table S5). Two IFN system genes (IFNphi3 and IFI44-like) and interleukin IL8 specifically raised their expression in clone F after CaHV infection, while CD antigen CD40 and chemokine CCL-c11a decreased their expression. Among these DEUs specific to d-A⁺ vs c-A⁺, approximately 18.1% (274/1510) up-regulated DEUs and 19.8% (628/3171) down-regulated DEUs were mapped into KEGG immune-related pathways or disease-related pathways (Additional file 8: Table S5). Different from a lot of IFN system genes specifically up-regulated in resistant clone H, only one IFN system gene, IRF1 was identified from the DEUs uniquely upregulated in susceptible clone A⁺. A lot of apoptosis or death related genes, such as tumor suppressor protein p53, bcl-2-like, Mcl1b, natural killer cell receptor 2B4 (CD244)-like isoform X1, X2 and X3, tumor necrosis factor receptor superfamily member tnfrsf9, death-associated protein kinase dapk3 and hematopoietic death receptor FAS, specifically increased their expression after CaHV infection. The results suggest that the susceptible clone A⁺ failed to protect against virus and evidently induced apoptosis or death. Additionally, interleukins (IL10 and IL22), interleukin receptors (IL6Rβ-like and IL21R), chemokines (CXC11, CCL-C25a and CCL-C24j) and CD antigens (CD27 and CD48) specifically raised their expression. In addition, several IFN system genes (IRF9 and IFITM3), complements (B/C2-A2, B/C2-A3, C3-H1, C3-S, C9 and CHFL4), chemokines (CXCb2, CXCL-C1 and CCL4), C-C chemokine receptors (CCRL2, CXCR1 and CXCR3), CD antigens (CD2, CD8beta-1a and CD86) and interleukin receptors (Il1r1, Il3r2 and IL2Rγ) were identified from the down-regulated DEUs specific to d-A⁺ vs c-A⁺.

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.

Fig. 7

Heatmap of IFN system genes. Hierarchical clustering is calculated by log₂ fold change values and the main clusters are lined out on the right. Up-regulated DEUs are shown in red and down-regulated DEUs in green

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.

In the battle between host and virus, the ability of host to mount a strong immune response is crucial for host survival [83, 84]. In fish, innate protective mechanisms include IFN, complement and cell-mediated cytotoxic activities [85]. As the first defense line of host, IFN system plays a vital role in resistance to viral infection [47]. From UV-inactivated GCHV-infected CAB cells, about 20 IFN system genes, such as TLR3, IFN, IRF1 [49], IRF3 [50], IRF7 [51], IRF9 [52], signal transducer and activator of transcription 1 (STAT1) [53], mediator of IRF3 activation (MITA) [54], Mx1 and Mx2 [55], PKR, PKZ [56], PKR-like [57], IFI58 and IFI56 [58], Gig1 and Gig2 [59], ISG15–1 and ISG15–2 [60], viperin [61], ubiquitin-specific protease 18 (USP18) [62] were identified. Similar to mammalian [86], fish initiates IFN response through recognizing viral products via TLRs and retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs), then triggers IRF3/7-dependent IFN response [10, 47, 87,88,89]. Interestingly, unique mechanisms underlying IRF3/7-dependent fish IFN response and auto-regulation of fish IFN gene expression were observed in UV-inactivated GCHV-infected CAB cells [47, 49, 50]. After exposing to CaHV, several pathways, “cytokine-cytokine receptor interaction” including IFNs, interleukins and chemokines, “Toll-like receptor signaling pathway” and “Jak-STAT signaling pathway” responsible for recognizing viral production or activating ISGs, were up-regulated in all three gibel carp clones (Fig. 4b; Additional file 8: Table S5), which clearly indicates that IFN response involves in common immune reaction of gibel carp to CaHV invasion.

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 IFN-mediated 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 head-kidney 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 tissue-specific 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 up-regulation 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 down-regulate 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 down-regulated 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 down-regulated 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 down-regulated 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.

Methods

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⁸ 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⁶ 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 real-time 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].