Transcriptome differences between susceptible and resistant fish
RNA sequencing was done on spleen samples of susceptible (S) and resistant (R) fish, each at two time points: pre-infection (day 0), and four days post infection (day 4). All three replicates of each treatment were highly correlated among themselves and thus, in a clustering analysis, based on normalized read counts of differentially expressed genes (DEGs), grouped together attesting for the high reproducibility of the data. Furthermore, results of the same fish type (susceptible or resistant) were more highly correlated (red color) to each other than pre- or post-infection results were (Fig 1A), demonstrating the distinct transcriptomes of these fish types. In the four comparisons between pairs of treatments (R0/S0, S4/R4, S4/S0 and R4/S4) a total of 2,369 DEGs were found. Relatively little overlap in DEGs between the four comparisons indicated that transcriptomes of the four treatments were considerably distinct (Fig 1B). For instance, only three genes were shared among all four comparisons, while 1,704 (72%) were comparison-specific. The largest overlap (435 DEGs) was found between S0/R0 and S4/R4, which represent fish-type differences regardless of the disease infection. The overlap between fish-types in response to infection (S4/S0 ∩ R4/R0) included only 69 DEGs.
Each of the four comparisons identified both up- and down-regulated genes relative to a reference treatment (the denominator of each comparison), all with a wide distribution of normalized read counts (Fig 2A). Some genes (n = 737) were more highly expressed in susceptible fish (Fig 2B), whereas other genes (n = 990) in resistant ones (Fig 2C). Only about ¼ of these genes (202/737 for susceptible and 258/990 for resistant) were shared between days 0 and 4, indicating that the transcriptional differences between fish types changed considerably in response to infection. These changes in expression between days 0 and 4, which reflect the response to infection, differed considerably between susceptible and resistant fish in two aspects. First, many more genes were differentially expressed in susceptible (n = 686) than in resistant (n = 199) fish. This broader response of susceptible fish applied to both up- and down-regulated genes (Fig 2D and E). Second, there was little overlap in the set of DEGs between susceptible and resistant fish. Only 9.3% (n=55) and 6.3% (n=14) of the up- and down-regulated DEGs, respectively, were shared between these fish types (Fig 2D and E). Therefore, in addition to baseline differences due to their different genetic background, these fish types had elicited a distinctive transcriptional response to CyHV-3 infection.
GO terms enrichment of DEGs in response to infection
The response to CyHV-3 infection is represented by the transcriptional differences elicited after infection and thus, the gene ontology (GO) of up-regulated gene lists between days 0 and 4 were studied (Fig 2D). In the 55 genes shared between susceptible and resistant fish the only significantly enriched GO term was ‘response to stress’ (GO:0006950), which reflects the general response of the fish to the challenge conditions. We further analyzed the 68 and 471 DEGs specific to the response of resistant and susceptible fish, respectively. The GO terms enriched in up-regulated DEGs were mostly immunity-related (Fig 3A). For susceptible fish, the most enriched term was ‘response to virus’ (GO:0009615) with 20 DEGs in carp that correspond to nine genes out of 48 included in this term based on the ZF annotation. The majority of these DEGs were interferons (IFNs) and interferon stimulated genes (ISGs). Another highly enriched term for susceptible fish was ‘chemokine-mediated signaling pathway’ (GO:0070098), with 11 DEGs in carp that correspond to seven out of 51 included in this term based on the ZF annotation (see Table 1 for a gene list of susceptible fish DEGs included in these terms). Other less enriched terms, including some not directly related to immune response, were also found for susceptible fish (Fig 3A).
Much fewer up-regulated DEGs were found in resistant fish, resulting in only three significantly enriched GO terms, all of which were immune response related (Fig 3A). These three terms were enriched also for susceptible fish, but, since tests were done on non-overlapping gene lists, these common terms were found to be enriched due to different genes, hence also their fold enrichment and test significance were different. For resistant fish, the highest fold enrichment was for ‘leukocyte migration’ (GO:0050900) with five carp genes that correspond to four out of 90 included in this term based on the ZF annotation. For susceptible fish, a lower fold enrichment of this term was found based on 12 carp genes that overlapped with those in the ‘chemokine-mediated signaling pathway’ term (GO:0070098).
Because ‘leukocyte migration’ was enriched in both fish types, we analyzed in more detail all DEGs in this term, including five DEGs specific to resistant, twelve to susceptible, and six shared by both (see Table 2 for a detailed gene list). Sixteen out of these 23 DEGs were chemokine genes, including six CC and ten CXC types. How many and which chemokines were up-regulated, differed between susceptible and resistant fish. Of the six CC chemokines, the five that were up-regulated in susceptible fish were ccl19 homologs and the one in resistant fish was a different chemokine, ccl35 (Fig 3B). Of the ten CXC chemokines, two homologs of cxcl11, two of cxcl18, one of cxcl20 and one of cxcl8, also known as interleukin-8 (il8), were up-regulated in susceptible fish, whereas in resistant fish, three homologs of cxcl8/il8 were up-regulated in response to infection. In addition, one homolog of cxcl20 was up-regulated in both fish types (Fig 3C). Of the seven non-chemokine genes in this GO term, five were up-regulated in both susceptible and resistant fish (Fig 3D). Altogether, these data provide further details on the general as well as specific immune response, differentiating resistant from susceptible fish. While susceptible fish mounted a broad response with clear signals of a viral infection, the response of resistant fish was associated with a specific chemokine-mediated, white blood cells migration response.
Expression pattern of duplicated chemokine genes
Chemokines comprise a large gene family with several homologs for many of the different family members in genomes of bony-fishes. Since the common carp is a tetraploid species, often times two paralogous copies are found in its genome per every ortholog in the genome of its closely related species, the diploid ZF. Since carp genes are annotated mostly after their homologs in ZF, in the case of gene families such as chemokines, separate carp genes (cypCars) are annotated with the same gene name from ZF. To better understand the possible contribution of chemokine homologs and of carp-specific paralogous gene pairs to resistance, the phylogenetic relationships between cypCar chemokines were studied. For all chemokines that one or more of their copies were found to be DEGs, DNA sequences of their transcripts from carp were extracted along with their orthologs in ZF and in red-bellied piranha as an outgroup. Sequences were aligned to construct a phylogenetic tree (Fig 4A). Based on how the evolutionary relationships among cypCars were resolved by the phylogenetic tree, differences in expression levels between carp copies within a gene clade were studied at three levels.
First, the proportion of DEGs out of the total gene copies was examined. Interestingly, different chemokine types had different proportions of their homologs changing expression in response to infection. For instance, 5/6 ccl19 homologs, 2/2 cxcl18 and 2/2 cxcl20 were DEGs. Conversely, only 1/4 ccl35 and 2/7 cxcl11 homologs were DEGs. Among the largest family, cxcl8/il8, 4/10 homologs were DEGs (Fig 4A).
Secondly, the specific expression differences of all ccl19 and cxcl8/il8 homologs in carp were tested using their normalized read counts. Significant differences in expression levels between genes within the ccl19 and cxcl8/il8 gene clades were found (all pairs Tukey-Kramer HSD, P < 0.05, Fig 4B). Among ccl19 homologs, four had a similar and very low expression level compared to two other homologs with significantly higher expression. Among cxcl8/il8 homologs, expression levels were diverse, ranging from nearly zero to almost 1,000 normalized read count. Thus, different chemokines diversified to a different extent in the expression of their homologs.
In summary, using a phylogentic approach to guide the transcriptional analysis revealed expression differences in response to infection among homologs of certain chemokines and even some functional divergence between carp-specific paralogs. Some of these transcriptional differences in homologs and paralogs were DEGs between susceptible and resistant fish.
Comparing the temporal response to infection
Since RNA-Seq was applied to susceptible and resistant fish only at day 4 post infection, it was interesting to study if susceptible and resistant fish elicited different responses or maybe similar responses but with different kinetics. Selected genes were chosen for RT-qPCR analysis in additional time points (day 2 and 8 post infection). Representing the ‘response to virus’ (GO:0009615) GO term, which was the most enriched one in susceptible fish, cypCar_00024055 [Radical S-adenosyl methionine domain containing 2 (rsad2) also known as Viperin-2 (vip2)] and cypCar_00017679 [Myxovirus resistance E (mxe)] were chosen. In the RNA-Seq analysis, both genes were up-regulated DEGs at day 4 in susceptible, but not in resistant fish. In the qRT-PCR analysis, no induction was observed at day 2, whereas at day 8 the genes were equally induced in both susceptible and resistant fish (Fig 5A and B).
Representing the GO term ‘leukocyte migration’ (GO:0050900), which was the most enriched one is resistant fish and was also enriched in susceptible fish, one ccl19a.1 copy (cypCar_19885) and two cxcl8/il8 copies (cypCar_47570 and cypCar_45160) were chosen. In the RNA-Seq analysis, cypCar_19885 (ccl19a.1) and cypCar_47570 (cxcl8/il8) were up-regulated DEGs at day 4 in susceptible, but not in resistant fish. In the qRT-PCR analysis, no induction was observed at day 2, whereas at day 8 the genes were equally induced in both susceptible and resistant fish (Fig 5C and D). However, cypCar_45160 (cxcl8/il8) that was up-regulated at day 4 in resistant, but not in susceptible fish, was not induced at all at day 2, and at day 8 it was kept induced in resistant and uninduced in susceptible fish (Fig 5E).
Thus, extending the RNA-Seq results to more time-points by RT-qPCR analyses indicated that several genes that were up-regulated at day 4 only in susceptible fish were later (day 8) also up-regulated in resistant fish. In contrast, the gene that was up-regulated at day 4 only in resistant fish, was not induced in susceptible fish. Taken together, the results suggested both qualitative and temporal differences between fish types in response to infection.
DEGs within previously identified QTLs
In our previous study [35], two novel QTLs associated with survival were identified and gene lists within QTL boundaries were extracted. Identifying DEGs within these QTLs can aid in highlighting candidate genes affecting survival. DEGs from all four pairwise comparisons were considered, since some of the differences in resistance might be attributed to genes with a basic difference in expression (day 0), while others to change in expression in response to the infection. From a few hundreds of genes located within the boundaries of each QTL1 and QTL2 intervals, 52 and 53 genes, respectively, were DEGs (Additional file 1: Table S1). Among these DEGs, some were annotated as immune related genes and were marked in the list as more promising candidate survival genes. Inside QTL1, a ccl20 homolog (cypCar_00050111) is more highly expressed in susceptible fish than in resistant fish at day 4. Additionally, inside QTL1, a cluster of DEGs annotated as fish specific tripartite motif-containing (TRIM) genes (finTRIMs – ftr) are present. Inside QTL2, a vitronectin b (vtnb) homolog (cypCar_00038793) is more highly expressed in resistant fish than in susceptible fish on day 0 and a Toll-like receptor-22 (tlr22) homolog (cypCar_00023948) is up-regulated between day 0 and 4 in resistant fish. Outside these QTLs, the previously identified candidate gene for survival, Interleukin-10a (il10a) [33, 35], was more highly expressed in susceptible fish than in resistant fish at day 4 (Additional file 1: Table S1). The expression of Interleukin-10b (il10b), the paralog of il10a, was similarly induced but the difference from resistant fish was insignificant. Thus, combining data of QTL mapping with differential expression and functional annotation highlighted candidate genes affecting disease survival.