Differences in responses of grass carp to different types of grass carp reovirus (GCRV) and the mechanism of hemorrhage revealed by transcriptome sequencing

Grass carp is an important farmed fish in China that is affected by serious disease, especially hemorrhagic disease caused by grass carp reovirus (GCRV). The mechanism underlying the hemorrhagic symptoms in infected fish remains to be elucidated. Although GCRV can be divided into three distinct subtypes, differences in the pathogenesis and host immune responses to the different subtypes are still unclear. The aim of this study was to provide a comprehensive insight into the grass carp response to different GCRV subtypes and to elucidate the mechanism underlying the hemorrhagic symptoms. Following infection of grass carp, GCRV-I was associated with a long latent period and low mortality (42.5%), while GCRV-II was associated with a short latent period and high mortality (81.4%). The relative copy number of GCRV-I remained consistent or decreased slightly throughout the first 7 days post-infection, whereas a marked increase in GCRV-II high copy number was detected at 5 days post-infection. Transcriptome sequencing revealed 211 differentially expressed genes (DEGs) in Group I (66 up-regulated, 145 down-regulated) and 670 (386 up-regulated, 284 down-regulated) in Group II. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed significant enrichment in the terms or pathways involved in immune responses and correlating with blood or platelets. Most of the DEGs in Group I were also present in Group II, although the expression profiles differed, with most DEGs showing mild changes in Group I, while marked changes were observed in Group II, especially the interferon-related genes. Many of the genes involved in the complement pathway and coagulation cascades were significantly up-regulated at 7 days post-infection in Group II, suggesting activation of these pathways. GCRV-I is associated with low virulence and a long latent period prior to the induction of a mild host immune response, whereas GCRV-II is associated with high virulence, a short latent period and stimulates a strong and extensive host immune response. The complement and coagulation pathways are significantly activated at 7 days post-infection, leading to the endothelial cell and blood cell damage that result in hemorrhagic symptoms.


Background
The grass carp (Ctenopharyngodon idellus) has been an important aquaculture species in China for over 60 years, accounting for more than 18% of total freshwater aquaculture production. The production of grass carp reached 5.5 million tons in 2014, making it the most highly consumed freshwater fish worldwide [1].
Grass carp hemorrhagic disease, caused by the grass carp reovirus (GCRV), is one of the most serious of these diseases [2]. GCRV, which was first isolated in China, belonging to the genus Aquareovirus of the family Reoviridae [3]. GCRV infects not only grass carp, but also rare minnow (Gobiocypris rarus), black carp (Mylopharyngodon piceus), and topmouth gudgeon (Pseudorasbora parva), causing hemorrhagic symptoms and death. Disease caused by GCRV outbreaks are frequent and result in huge economic losses in the aquaculture industry. Consequently, GCRV is of particular interest to fish breeding geneticists aiming to identify strategies for disease-resistant breeding [4][5][6][7][8][9].
Recently, the genome sequences of a number of GCRV strains isolated in China have been determined [10][11][12]. Sequence comparisons and analysis showed that GCRV could be divided into three distinct subtypes. GCRV-873, which is a representative strain of type I (GCRV-I), infected C. idellus kidney (CIK) cells and induced obvious cytopathic effects (CPE). GCRV-HZ08, which is a representative strain of type II GCRV (GCRV-II), resulted in 80% mortality in yearling fish, while no obvious CPEs were observed in GCRV-II infected CIK cells. Type III GCRV (GCRV-III) is not widely distributed in China and only one strain was found (GCRV-104), which also induced CPE in CIK cells.
Of the three types of GCRV, GCRV-873 was the first fish virus to be characterized and sequenced [13] and was used as the target in early studies focusing on disease-resistant breeding and virus-host interactions. However, recent studies showed that most of the GCRV isolated in Southern China are type II GCRVs, such as GCRV-HZ08, GCRV-GD108 and GCReV-109 [10][11][12]. GCRV types I and II both cause hemorrhagic disease in grass carp, although the hemorrhage mechanism is unknown. Moreover, these virus types have significantly different nucleotide sequences, viral-encoded protein structures, and pathogenicity in grass carp and CIK cells. Although many studies on GCRV have been conducted, most were restricted to investigation of the virus itself and differences in pathogenesis and immune responses to different types of GCRV in grass carp remain to be elucidated.
In this study, grass carp were infected with two types of GCRV (GCRV-I and GCRV-II) and the pathogenesis was investigated by transcriptome sequencing, real time quantitative PCR (RT-qPCR) and mortality rates. Details of the transcriptional events following GCRV infection have been reported previously [14]; therefore, the aim of the current study was to provide a comprehensive insight into the responses of grass carp to different types of GCRV and to reveal the mechanism underlying the hemorrhagic symptoms. Our study will provide guidance for the development of novel vaccines and diseaseresistant breeding of grass carp.

Viruses
Two types of grass carp reovirus were used in the study. One GCRV strain, isolated in Honghu City, Hubei Province, China in May 2015, was classified as a type I GCRV due to the high similarity (97.3%) of the S2 segment to the typical type I strain, GCRV-873. Another GCRV strain, isolated in Huanggan city, also in Hubei Province, in July 2015, was classified as a type II GCRV due to the high similarity (98.4%) of the S2 segment to the typical type II strain, GCRV-HZ08. The two types of GCRV were designated GCRV-I and GCRV-II for the purposes of this study and both were diluted to the same titer (2.97 × 10 3 RNA copy/μl) for use in experiments.

Experimental fish
Healthy full-sib grass carp were used in the study at 3 months of age, weighing 3-5 g and with an average length of 8 cm. The fish were obtained from the Guan Qiao Experimental Station, Institute of Hydrobiology, Chinese Academy of Sciences, and acclimatized in aerated fresh water at 26-28°C for one week before processing. Fish were fed with a commercial diet twice a day and water was exchanged daily. If no abnormal symptoms were observed, grass carp were selected for further study. Fish were then divided into three groups (approximately 150 per group) that were maintained in separate tanks.

Virus challenge experiment and sample collection
After no abnormal symptoms were observed in the three groups, virus challenge experiments were carried out. Fish in the Groups I and II were infected with 200 μl GCRV-I (2.97 × 10 3 RNA copy/μl) or 200 μl GCRV-II (2.97 × 10 3 RNA copy/μl) by intraperitoneal injection, respectively, while fish from Group III were injected with 200 μl PBS as a Control group. At 1, 3, 5, and 7 days post-injection, 15 fish that contained three biological duplicates (n = 5 for each biological duplicate) from each group were collected and the kidneys were removed for analysis . The samples were designated I-1, I-3, I-5, I-7,  II-1, II-3, II-5, II-7, c-1, c-3, c-5, and c-7 (three biological duplicates for each sample). The remaining fish were monitored carefully and the number of dead fish in each group was counted every day. The experiment was concluded and the total mortality was calculated when no mortality was recorded for seven consecutive days.

RNA isolation, library construction and sequencing
RNA of kidneys that collected above was isolated using TRizol reagent (Invitrogen, USA) according the manufacturer's protocol. RNA concentration was measured using the Qubit RNA assay kit (Life Technologies, USA), and integrity was assessed with the RNA Nano 6000 assay kit (Agilent Technologies, USA). RNA of sufficiently high quality was used in library construction. Sequencing libraries were generated using the NEBNext Ultra RNA library prep kit for Illumina (New England Biolabs, USA) following the manufacturer's protocol. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads and fragmented by the NEBNext first strand synthesis reaction buffer (New England Biolabs). First strand cDNA was synthesized using a random hexamer primer and M-MuLV reverse transcriptase. Second strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. After adenylation of the 3′ end of DNA fragments, NEBNext adaptors with a hairpin loop structure were ligated in preparation for hybridization. Subsequently, 3 μL USER enzyme (New England Biolabs, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C prior to PCR using Phusion High-fidelity DNA polymerase, universal PCR primers and index (X) primer. Finally, PCR products were purified using an AMPure XP system and library quality was assessed using an Agilent Bioanalyzer 2100 system. Libraries were sequenced on an Illumina Hiseq X Ten platform and 150 bp pair-end reads were generated.

Data analysis
Raw data reads in fastq format were initially processed using in-house perl scripts. In this step, clean data (clean reads) were obtained by removing adapter, poly-N and poor quality data. The Q20, Q30, and GC contents of the clean data were calculated, and all downstream analysis was performed using the clean high quality data.
Clean data were mapped to the grass carp reference genome [15] using TopHat2 software [16]. Allowing for two base mismatches in the mapping process, total mapped reads were calculated, and the mapped regions (exon, intron, and intergenic) were counted.
HTSeq software was used to count the number of reads mapped to each gene [17] and the reads per kilobase of the exon model per million mapped reads (RPKM) were calculated for each gene based on the length of the gene and the number of reads mapped to the gene [18].

Differential expression analysis
Differential expression analysis of two groups/conditions was performed using the DESeq package [19]. The resulting P-values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05 (q value <0.05) in DESeq analysis were assigned as differentially expressed genes (DEGs).
Gene Ontology (GO) annotation of the genes was performed using ClueGO and CluePedia [20,21]. In GO enrichment analysis, only categories with a low P-value (P < 0.05) were considered as enriched in the network as determined by two-sided hypergeometric statistical tests employing the Benjamini and Hochberg approach to false discovery rate correction.
The Kyoto Encyclopedia of Genes and Genomes (KEGG) database is used to provide high-level functional information on biological systems of molecules, cells, organisms and ecosystems, and is particularly powerful for the evaluation of large-scale molecular datasets generated by genome sequencing and other high-throughput experimental approaches [22]. In this study, KOBAS software was employed to test the statistical enrichment of DEGs in KEGG pathways [23]. KEGG terms with corrected P < 0.05 were considered to indicate statistical significance.

Validation of DEGs by RT-qPCR
To confirm the reliability of data obtained by RNA-seq, 10 DEGs were selected for validation by RT-qPCR. The primers are listed in Additional file 1. Only primers with efficiency of 90%-110% were used for RT-qPCR analysis. First strand cDNAs were obtained using a random hexamer primer and the ReverTra Ace kit (Toyobo, Japan). RT-qPCR was carried out using a fluorescence quantitative PCR instrument (Bio-Rad, USA). Each RT-qPCR mixture contained 0.8 μL forward and reverse primers (for each primer), 1 μL template, 10 μL 2× SYBRgreen master mix (TOYOBO, Japan), and 7.4 μL ddH 2 O. Three replicates were included for each sample and the β-actin gene was used as an internal control to for normalization of gene expression. The program for RT-qPCR was as follows: 95°C for 10 s, 40 cycles of 95°C for 15 s, 55°C for 15 s, and 72°C for 30 s. Relative expression levels were calculated using the 2 -△△Ct method [24]. Data represent the mean ± standard deviation of three replicates.

Statistical analysis
The statistical significance between Group I and Group II was determined by one-way ANOVA. Differences were considered significant at P < 0.05. P < 0.05 was denoted by *.

Mortality of grass carp infected with GCRV-I or GCRV-II
The mortality curves of the three groups are shown in Fig. 1. In the GCRV-I-infected group, the total mortality of 42.5% was reached at 15 days post-infection, with the first fish dying at 10 days post-infection. However, in the GCRV-II-infected group, the total mortality of 81.4% was reached at 15 days, with the first death recorded as early as 8 days post-infection. In the Control group, two dead individuals were observed, giving a total mortality of 2.2%. Moreover, fish that died after infection with GCRV-I or GCRV-II showed hemorrhagic symptoms, especially in the muscle, whereas no hemorrhagic symptoms were observed in the Control group fish (Fig. 2).

Relative copy numbers of the two types of GCRV in grass carp
To determine dynamic changes in the levels of the two types of GCRV in infected fish, the relative copy numbers of the viruses were examined by RT-qPCR using specific primers for the S6 segments of the two types of GCRV. For convenience, the relative copy number of GCRV in one day post-infection in Group I was used as a reference for normalization. As shown in Fig. 3, the relative copy number of GCRV-I in Group I was consistent or decreased slightly throughout the experiment. In contrast, the relative copy number of GCRV-II in Group II was extremely low at 1 and 3 days post-infection, followed by a marked increase at 5 days post-stimulation and a reduction to levels similar to those of GCRV-I in Group I at 7 days post-infection. Obviously, the two types of GCRV showed different dynamic curves during infected fish.

Preliminary analysis of transcriptome sequencing data
To further reveal the mechanism underlying the differences in sensitivity of grass carp to the two types of GCRV, we performed RNA-seq analysis on samples collected from three groups at different time-points post-infection. Three duplicates of each sample were processed, yielding a total of 36 libraries, which were sequenced on an Illumina Hiseq X Ten platform to generate 150 bp pair-end reads. As shown in Table 1, the raw reads, clean reads, clean base Q20, Q30, and   (Table 1). All libraries gave Q20 ≥ 95%, Q30 ≥ 87%, and mapped percent ≥87%. These results confirmed the high quality of the sequencing data and suitability for further analysis. The sequencing data in this study have been deposited in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (accession number: SRP095827).

Identification of DEGs
DEGs among these samples were identified by subjecting the data to a series of paired-comparisons. In the analysis, samples from Group I (I-1, I-3, I-5, and II-7) and Group II (II-1, II-3, II-5, and II-7) were compared with samples that from the Control group at the corresponding time-points the same time. The numbers of DEGs identified from the different paired-comparisons are listed in Table 2. Comparisons with the Control group revealed 66 up-regulated and 145 down-regulated genes in Group I, whereas 386 up-regulated and 284 down-regulated genes were identified in Group II. Venn diagram analysis of the 49 upregulated and 115 down-regulated genes found in both groups is shown in Fig. 4. In detail, 25, 4, 38 and 16 genes were un-regulated, whereas 102, 13, 35, and 28 genes were down-regulated in Group I at 1, 3 5 and 7 days postinfection, respectively. In Group II, 42, 94, 203, and 139 genes were up-regulated and 29, 54, 234, and 29 genes were down-regulated at 1, 3 5 and 7 days post-infection, respectively. Detailed information of these DEGs is shown in Additional file 2. The DEGs that could not be functionally annotated are listed as "unknown".

GO and KEGG enrichment analysis of the DEGs
GO enrichment analysis was performed to investigate the possible roles of these DEGs. Of the three main categories (biological process, molecular function, and cellular component), most of the GO terms belonged to the biological process category, suggesting the occurrence of a series of molecular events in grass carp after GCRV infection. Many of the significantly enriched GO terms in both groups were involved in immune responses, such as defense response, response to organic substance, inflammatory response, acute inflammatory response, and innate immune response, indicating that grass carp respond strongly to both types of GCRV infection. Moreover, some GO terms associated with blood or platelets, such as blood microparticles, blood coagulation, platelet activation, and platelet degranulation, were also significantly enriched. The top five enriched GO terms in both groups are listed in Table 3 and details of the GO terms are shown in Additional file 3.
KEGG enrichment analysis was also performed for the DEGs in both groups. In general, many the significantly enriched KEGG terms were also involved in immune responses, such as complement and coagulation cascades, Staphylococcus aureus infection, the peroxisome proliferator-activated receptor (PPAR) signaling pathway, graft-versus-host disease, and allograft rejection. In addition, some terms involved in metabolism and biosynthesis were also significantly enriched. These terms included glycolysis/gluconeogenesis, phenylalanine, tyrosine and tryptophan biosynthesis, as well as vitamin digestion and absorption, and fat digestion and absorption. The top five enriched KEGG terms in both groups are also listed in Table 3 and details of the KEGG terms are shown in Additional file 4.

Expression patterns of DEGs in the key GO terms
The GO terms involved in blood coagulation, defense response, innate immune response, inflammatory response, apoptotic process, and metal ion transport were selected for further investigation. Comparisons at each time-point after infection revealed significant differences in the expression patterns of these DEGs in the two groups (Fig. 5). The DEGs in Group I showed mild change in expression at the four time-points, with log 2 fold changes ≤2 for most of genes. However, a significant variation in the expression levels of the DEGs in Group II was observed, with the most marked variation at 5 days post-infection. Many of DEGs in Group II showed log 2 fold changes ≥2 or even ≥4.

Identification of DEGs shared between the two groups
Venn diagram analysis showed that the two groups shared 49 up-regulated and 118 down-regulated DEGs, regardless of the time-points post-infection. These shared Fig. 3 Relative copy number of GCRV in Group I and Group II. The relative copy numbers of GCRV in Groups I and II were examined by using specific primers for the S6 segment. The copy number of GCRV at 1 day post-infection in Group I was used as a reference for normalization. Data represent mean ± standard deviation of three replicates. Significant difference (P < 0.05) between the two groups was indicated with asterisks (*) DEGs were subjected to further analysis to investigate the common events after infection with different types of GCRV. Of the 49 up-regulated DEGs, many were involved in anti-viral immune response, blood coagulation, response to stress, protein degradation, antigen presentation, and transcription/translation. Although these DEGs were up-regulated in both groups, the expression patterns differed. The representative upregulated DEGs are listed in Table 4. Of the down-regulated DEGs, many were involved in in blood coagulation, cytoskeleton, complement activation, iron transport, inflammatory response, and metabolism. The representative up-regulated DEGs are also listed in Table 4. Interestingly, the representative DEGs (both up-regulated and down-regulated) belonging to the same categories showed similar expression patterns in both groups. Details of these DEGs in both groups are shown in Additional file 5.  Despite the many differences in the up-regulated DEGs between the two groups, many were categorized as immune response, response to stress, antigen presentation, and blood coagulation. In addition, most of the DEGs in Group I showed mild upregulation, with log 2 fold changes ≤2, whereas those in Group II were intensely up-regulated, with log 2 fold changes ≥2. The DEGs with log 2 fold changes ≥2 in Group I and ≥3 in Group II are shown in Table 5. There were fewer downregulated DEGs in Group II than in Group I. Many of the DEGs were involved in iron transport, transcription/translation, cytoskeleton, enzyme activity, and binding activity, suggesting inhibition of the host transcription/ translation machinery by GCRV. Interestingly, many genes encoding the middle subunit of ferritin were significantly downregulated. The DEGs with log 2 fold changes ≤ − 2 in both groups are shown in Table 5.

DEGs involved in complement and coagulation cascades
KEGG enrichment analysis showed that the "complement and coagulation cascades" pathway was the most significantly enriched in both groups. The expression patterns of the DEGs involved in these pathways were further investigated to elucidate the mechanism of hemorrhage in infected fish (Fig. 7). Interestingly, most of the DEGs showed similar expression patterns in each group, with the exception of the urokinase plasminogen activator surface receptor (PLAUR) gene, which showed an opposing expression pattern compared with other DEGs. The DEGs in Group I showed only slight changes in expression, with it increased levels detected at 3 days and decreased expression at 1, 5, and 7 days postinfection. Most of the DEGs in Group II showed mild changes at 1 and 3 days post-infection. However, significantly decreased expression of the DEGs was observed at 5 days post-infection, followed by marked upregulation at 7 days post-stimulation, suggesting activation of the "complement and coagulation cascades" pathways.

Confirmation of DEGs by RT-qPCR
To confirm the RNA-seq data, a total of 10 DEGs (5 involved in "immune response" and 5 involved in "complement and coagulation cascades") were selected for RT-qPCR analysis. These genes included Mx-1, IRF3, IRF7, CCL7, IFIT1, PLAUR-1, alpha-2-macroglobulin 3 (A2M-3), complement C3c alpha chain fragment 1 (C3-1), coagulation factor IX (F9), and complement component 9 (C9). The relative expression levels of DEGs in Groups I and II were calculated as the ratio of gene expression levels relative to those in the Control group at the corresponding time-point. As shown in Fig. 8, overall, the expression patterns of all 10 DEGs identified by qPCR were similar to those obtained in RNA-seq analyses, although the relative expression levels were not completely consistent. Moreover, RT-qPCR also showed that the PLAUR-1 gene expression pattern was opposite to that of the other genes that involved in "complement and coagulation cascades". Therefore, the results of the RT-qPCR analysis confirmed the reliability and accuracy of the RNA-seq data.

Discussion
In the study, we used RNA-seq to elucidate the the mechanism of hemorrhage and the different immune responses induced in grass carp infected with GCRV-I and GCRV-II. Mortality analysis showed that grass carp infected with GCRV-I began to die at 10 days post-infection with a total mortality of 42.5%, whereas the first grass carp infected with GCRV-II died at 8 days postinfection and the total mortality was much higher, at 81.4%. The results suggested that the GCRV-II subgroup is associated with high virulence and a short latent period, while the GCRV-I subgroup is associated with low virulence and a longer latent period. Interestingly, changes in the relative copy numbers of the two types of GCRV during the course of infection showed results consistent with the theory mentioned in the above. The copy number of GCRV-I remained stable or even decreased slightly from day 1 to day 7 post-infection, whereas GCRV-II showed a rapid increase in copy number at 5 days post-stimulation. High virulence, short latent period, and high mortality could accelerate the spread of virus [30]. Our results explained that why the GCRV isolated recently in China belonged to GCRV-II rather than GCRV-I [10][11][12]. Moreover, the number of DEGs in the two groups is  Mx-1, IRF2, IRF3, IRF7, CCL7,  CCL19, IL8, IL11, IFIT1, and IFIH1 in both groups also consistent with the relative copy numbers. It can be speculated that fewer DEGs were identified in Group I due to the low copy numbers at the four timepoints, while numerous DEGs were discovered in Group II, especially at 5 days post-infection due to the high GCRV-II copy number of at this time-point. The consistency of these data suggests the reliability of our experimental design. Furthermore, the relatively low number of DEGs found in Group I compared with Group II suggest the induction of a strong and extensive response in Group II compared with the relatively mild response in Group I. The changes in DEGs in Group II were more dramatic than those in Group I. The variation in the expression pattern of DEGs among the selected GO terms was also greater in Group II than that in Group I. However, most of the DEGs in Group I were also identified in Group II, although the differences in expression patterns were consistent with the variation in the course of infection between the two groups. Therefore, we hypothesized that GCRV-I and GCRV-II induced similar responses in grass carp, although the response to GCRV-II was more robust and extensive. Previous studies revealed similar results in rainbow trout, which showed similar responses to high and low virulence strains of infectious hematopoietic necrosis virus (IHNV), although the level of the response and DEG profiles differed [31,32].
Most of the up-regulated DEGs common to Groups I and II were associated with functions in anti-viral immune responses, blood coagulation, response to stress, protein degradation, antigen presentation, and transcription/translation, suggesting that a series of events, especially immune responses, occur after infection regardless of the GCRV subtype. Many of the down-regulated DEGs common to both Groups I and II were found to be involved in cytoskeleton, metabolism, iron transport, enzyme activity, and binding activity. This association was particularly marked in Group I at 1 day postinfection. The association of up-regulated DEGs with protein degradation, while down-regulated DEGs were associated with metabolism, enzyme activity, and binding activity suggests that the host translation machinery is hjjacked or shut-down by GCRV to facilitate the replication and spread of the virus. A similar phenomenon has been observed in other fish after virus infection [33]. Interestingly, DEGs involved in cytoskeleton and iron transport were down-regulated in both Groups I and II. These observations are consistent with reports that actin cytoskeletal dynamics are important for T lymphocyte activation and migration [34,35] and the important role of iron homeostasis in the host response to infection [36][37][38].
GCRV causes hemorrhagic symptoms in infected fish, although the mechanism remains unknown. Interestingly, Fig. 7 Heatmap of DEGs involved in the complement pathway and coagulation cascades. Heatmap showing log 2 fold changes in differentially expressed genes in both groups involved in the complement pathway and coagulation cascades GO enrichment analysis revealed that some GO terms associated with blood or platelets, such as blood microparticles, blood coagulation, platelet activation, and platelet degranulation were significantly enriched, especially in Group II, indicating the existence of a correlation between these GO terms and the hemorrhagic symptoms. Moreover, KEGG analysis showed that DEGs associated with "complement and coagulation cascades" were the most Fig. 8 Confirmation of the transcriptome sequencing data by RT-qPCR. Five DEGs participating in immune responses and five DEGs involved in "complement and coagulation cascades" were selected for RT-qPCR. The relative expression levels of DEGs was calculated as the ratio of gene expression level relative to that in the Control group at the same time-point. All data represent the mean ± standard deviation of three replicates