Transcriptomics Sequencing Provides Insights into Understanding the Mechanism of Grass Carp Reovirus Infection

Grass carp is an important aquaculture fish species in China that is affected by severe diseases, especially haemorrhagic disease caused by grass carp reovirus (GCRV). However, the mechanisms of GCRV invasion and infection remain to be elucidated. In the present study, Ctenopharyngodon idellus kidney (CIK) cells were infected with GCRV, harvested at 0, 8, 24, and 72 h post infection, respectively, and then subjected to transcriptomics sequencing. Each sample yielded more than 6 Gb of clean data and 40 million clean reads. To better understand GCRV infection, the process was divided into three phases: the early (0–8 h post infection), middle (8–24 h post infection), and late (24–72 h) stages of infection. A total of 76 (35 up-regulated, 41 down-regulated), 553 (463 up-regulated, 90 down-regulated), and 284 (150 up-regulated, 134 down-regulated) differently expressed genes (DEGs) were identified during the early, middle, and late stages of infection, respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed that DEGs were mainly involved in carbohydrate biosynthesis, transport, and endocytosis in the early stage, phagocytosis and lysosome pathways were mainly enriched in the middle stage, and programmed cell death, apoptosis, and inflammation were largely associated with the late stage. These results suggest GCRV infection is a gradual process involving adsorption on the cell surface, followed by endocytosis into cells, transport by lysosomes, and eventually resulted in cell necrosis and/or apoptosis. Our findings provide insight into the mechanisms of grass carp reovirus infection.


Introduction
The grass carp (Ctenopharyngodon idellus) has been reared in China for more than 60 years as an important aquaculture species. The production of grass carp reached 5.8 million tons in 2015, accounting for more than 13% of the world's freshwater aquaculture production [1]. Grass carp haemorrhage disease, caused by grass carp reovirus (GCRV), is one of the most damaging diseases, resulting in huge economic losses to the aquaculture industry of grass carp [2]. GCRV was first isolated in China and belongs to the genus Aquareovirus of the family Reoviridae [3]. GCRV not only infects grass carp but also infects the rare minnow (Gobiocypris rarus), black carp (Mylopharyngodon piceus) and topmouth gudgeon (Pseudorasbora parva), causing haemorrhagic symptoms and death in these species.
Grass carp haemorrhage disease outbreaks are frequent and result in huge economic losses in the aquaculture industry, hence the GCRV is of particular interest to geneticists aiming to develop strategies Grass carp haemorrhage disease outbreaks are frequent and result in huge economic losses in the aquaculture industry, hence the GCRV is of particular interest to geneticists aiming to develop strategies for breeding disease-resistant fish. Recent studies on GCRV have mainly focused on identification of genes involved in immune responses to GCRV infection [4,5], cloning and characterization of GCRV-encoded proteins [6,7], the development of vaccines against GCRV [8,9], and transcriptomic sequencing of discovered immune-related genes or uncovered the antiviral immunity mechanism of grass carp infected with GCRV [10][11][12][13][14]. However, the mechanisms of grass carp reovirus invasion and infection remain poorly explored.
In the present study, C. idellus kidney (CIK) cells were infected with GCRV, harvested at 0, 8, 24, and 72 h post infection, respectively, and subjected to transcriptomics sequencing. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed on differentially expressed genes (DEGs) identified during different stages of infection. Our findings expand our understanding of the mechanisms underlying grass carp reovirus infection.

GCRV Infection Induces Cytopathic Effects
As shown in Figure 1, no cytopathic effect (CPE) was observed in mock-infected cells. However, in cells infected with GCRV, a CPE was observed as early as 8 h post infection, and it became more pronounced with over time. At 72 h post infection, it could be seen that most cells were necrotic.

Cell Counting Kit-8 Assay
As shown in Figure 2, the mock-infected cells showed the viability near 100%. However, for cells that infected for 8, 24, and 72 h, only 74.66%, 56.88%, and 36.95% of cell viability were observed. Therefore, the result suggested that GCRV infection is efficiency.

Cell Counting Kit-8 Assay
As shown in Figure 2, the mock-infected cells showed the viability near 100%. However, for cells that infected for 8, 24, and 72 h, only 74.66%, 56.88%, and 36.95% of cell viability were observed. Therefore, the result suggested that GCRV infection is efficiency. Figure 2. Cell viability assay after infected with GCRV. After 0, 8, 24 and 72 h post infection, cells were evaluated using CCK-8 assay, respectively. The mean of three replicates was shown with the ±standard deviation (S.D.). Significant difference between the control and treated group was indicated with asterisks (*: p < 0.05; **: p < 0.01)

Preliminary Analysis of Transcriptomics Sequencing Data
As shown in Table 1, the raw reads, clean reads, clean bases, Q20, Q30, and mapped percentage were recorded for each library. For all libraries, Q20 ≥ 96%, Q30 ≥ 92% and mapped percentage ≥87%. These results showed the high quality of the sequencing data and ensured suitability for further analysis. All sequencing data have been uploaded to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI; accession number SRP111072). Moreover, principal component analysis (PCA) was performed on each sample according to the expression level. As shown in Figure 3, mock-infected cells were clearly distinct from cells infected with GCRV. For samples infected with GCRV, correlation values were proportional to the time post-infection; samples from 8 h did not cluster with those from 24 and 72 h. Therefore, the results suggested that infection efficiency and differences existed in each sample post infection.  were evaluated using CCK-8 assay, respectively. The mean of three replicates was shown with the ±standard deviation (S.D.). Significant difference between the control and treated group was indicated with asterisks (*: p < 0.05; **: p < 0.01)

Preliminary Analysis of Transcriptomics Sequencing Data
As shown in Table 1, the raw reads, clean reads, clean bases, Q20, Q30, and mapped percentage were recorded for each library. For all libraries, Q20 ≥ 96%, Q30 ≥ 92% and mapped percentage ≥87%. These results showed the high quality of the sequencing data and ensured suitability for further analysis. All sequencing data have been uploaded to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI; accession number SRP111072). Moreover, principal component analysis (PCA) was performed on each sample according to the expression level. As shown in Figure 3, mock-infected cells were clearly distinct from cells infected with GCRV. For samples infected with GCRV, correlation values were proportional to the time post-infection; samples from 8 h did not cluster with those from 24 and 72 h. Therefore, the results suggested that infection efficiency and differences existed in each sample post infection.

Identification of Differently Expressed Genes (DEGs)
To better understand GCRV infection, the process was divided into three phases: early (0-8 h post infection), middle (8-24 h post infection), and late (24-72 h) stages. In each phase, data from the latter time point were compared with the former (8 vs. 0, 24 vs. 8, and 72 vs. 24) to identify differentially expressed genes (DEGs). As shown in Table 2, a total of 76 (35 up-regulated, 41 downregulated), 553 (463 up-regulated, 90 down-regulated), and 284 (150 up-regulated, 134 downregulated) DEGs were identified during early, middle, and late stages of infection, respectively. Detailed information related to these DEGs is shown in Table S1. DEGs that could not be functionally annotated are listed as 'unknown'.

GO and KEGG Enrichment Analysis
Go enrichment analysis was performed to investigate the putative roles of DEGs. The early stage of infection was mainly related to iron transport, carbohydrate metabolism, and inclusion bodies. For the middle stage of infection, iron ion binding, lysosome categories were significantly enriched. Meanwhile, the late stage was mainly associated with metabolic processes such as lipid biosynthetic,

Identification of Differently Expressed Genes (DEGs)
To better understand GCRV infection, the process was divided into three phases: early (0-8 h post infection), middle (8-24 h post infection), and late (24-72 h) stages. In each phase, data from the latter time point were compared with the former (8 vs. 0, 24 vs. 8, and 72 vs. 24) to identify differentially expressed genes (DEGs). As shown in Table 2, a total of 76 (35 up-regulated, 41 down-regulated), 553 (463 up-regulated, 90 down-regulated), and 284 (150 up-regulated, 134 down-regulated) DEGs were identified during early, middle, and late stages of infection, respectively. Detailed information related to these DEGs is shown in Table S1. DEGs that could not be functionally annotated are listed as 'unknown'.

GO and KEGG Enrichment Analysis
Go enrichment analysis was performed to investigate the putative roles of DEGs. The early stage of infection was mainly related to iron transport, carbohydrate metabolism, and inclusion bodies. For the middle stage of infection, iron ion binding, lysosome categories were significantly enriched. Meanwhile, the late stage was mainly associated with metabolic processes such as lipid biosynthetic, steroid metabolic, and alcohol biosynthetic metabolic. The top five enriched GO terms in each of the three stages are listed in Table 3, and details of the GO terms are included in Table S2. KEGG enrichment analysis was also performed for DEGs identified in the three stages. During the early stage of infection, DEGs were mainly involved in pathways related to focal adhesion, ECM (extracellular matrix)-receptor interactions, cytokine-cytokine receptor interactions. In the middle stage, pathways related to virus entry into cells and cell recognition were enriched, including the phagosome, lysosome, and hippo signaling pathways. In the late stage of infection, the most enriched pathways were the biosynthesis of steroids and the terpenoid backbone, and degradation of ketone bodies. The top five enriched KEGG terms in each of the three stages are listed in Table 4, and details of the KEGG terms are included in Table S3.

Expression of DEGs in Key KEGG Pathways
Expression patterns of DEGs in the key KEGG pathways ECM-receptor interactions, focal adhesion, phagosome, and lysosome KEGG pathways were further evaluated. As shown in Table 5, DEGs related to ECM-receptor interactions and focal adhesion, such as col5a1, itgb7, and sptbn5 was up-regulated during the early stage of infection but was down-regulated during the middle and late stages. By contrast, expression of DEGs involved in phagosomes and lysosomes, such as itgb2, srb1, and sel1l, followed the opposite trend, and was down-regulated during the early stage of infection, but up-regulated during the middle and late stages. Comparison of the four KEGG pathways at each time point after infection revealed a significant difference in the expression patterns of these DEGs ( Figure 4). Table 5. DEGs associated with four major pathways.

Pathway
Gene Name Log2 Fold Change

Key DEGs Identified in the Three Stages of Infection
The more significant DEGs may play important roles in the responses to changes in the environment, hence they were identified and annotated. As shown in Table 6, most of the upregulated genes between 0 and 8 h post infection such as protein tyrosine phosphatase epsilon (ptpε), integrin alpha 7 (itga7) and periplakin (ppl), are associated with cell signaling, whereas downregulated genes such as C-X-C motif chemokine 10-like (cxcl10) and interferon alpha-inducible protein 27-like (ifi27l) are associated with immunity. Genes up-regulated between 8 and 24 hpi are primarily immune-related, including viperin, antigen peptide transporter 1 (tap1), toll-like receptor 3 (tlr3), and interferon regulatory factor-1 (irf11). Meanwhile, down-regulated genes such as angiopoietin-related protein 2-like (angptl2), phosphoglycolate phosphatase (pgp), and RAP1 GTPase activating protein 1 (rap1gap1) are mainly related to blood vessels, which may explain why haemorrhage symptoms were observed in grass carp following GCRV infection. Other downregulated genes include protein damage repair-related genes such as methionine sulfoxide oxidase 3b (mical3b) and heat shock protein beta 11 (hspb11). Finally, between 24 and 72 h, up-regulated genes are mainly related to programmed cell death, apoptosis, and inflammation, consistent with the high level of cell viability reduction observed at 24 and 72 hpi (Figure 2). The main down-regulated genes

Key DEGs Identified in the Three Stages of Infection
The more significant DEGs may play important roles in the responses to changes in the environment, hence they were identified and annotated. As shown in Table 6, most of the up-regulated genes between 0 and 8 h post infection such as protein tyrosine phosphatase epsilon (ptpε), integrin alpha 7 (itga7) and periplakin (ppl), are associated with cell signaling, whereas down-regulated genes such as C-X-C motif chemokine 10-like (cxcl10) and interferon alpha-inducible protein 27-like (ifi27l) are associated with immunity. Genes up-regulated between 8 and 24 hpi are primarily immune-related, including viperin, antigen peptide transporter 1 (tap1), toll-like receptor 3 (tlr3), and interferon regulatory factor-1 (irf11). Meanwhile, down-regulated genes such as angiopoietin-related protein 2-like (angptl2), phosphoglycolate phosphatase (pgp), and RAP1 GTPase activating protein 1 (rap1gap1) are mainly related to blood vessels, which may explain why haemorrhage symptoms were observed in grass carp following GCRV infection. Other down-regulated genes include protein damage repair-related genes such as methionine sulfoxide oxidase 3b (mical3b) and heat shock protein beta 11 (hspb11). Finally, between 24 and 72 h, up-regulated genes are mainly related to programmed cell death, apoptosis, and inflammation, consistent with the high level of cell viability reduction observed at 24 and 72 hpi (Figure 2). The main down-regulated genes are involved in steroid synthesis, while others are associated with microtubules and blood vessels. Details of these DEGs in the three stages are included in Table S4.

Expression Patterns of Key Interferon-Related Genes
The expression patterns of key interferon-related genes were examined in this study. As shown in Figure 5, except gene mx-1, most of the genes showed log2 fold change ≥2 or even ≥10 at 24 or 72 h post-infection, whereas the expression level decreased at 8 h post-infection. These results suggested that strong immune response induced by GCRV infection.

Discussion
At present, previous studies on grass carp transcriptome profiles in response to grass carp reovirus infection have been carried out [10][11][12][13][14]. Most of them involved in the discovery of immune-related genes or uncovered the antiviral immunity mechanism of grass carp. However, the purpose of our study is to understand the whole process of grass carp reovirus infection.
In this study, we used the CIK cells as an in vitro model to explore the molecular changes following reovirus infection. CIK cells were mock-infected or infected with GCRV and harvested at 0 (mock), 8,24, and 72 hpi. To better understand the process of GCRV entry into cells, it was divided into three stages: early (0-8 h), middle (8-24 h), and late (24-72 h). In each stage, data from the latter time point were compared with that from the former to identify DEGs. The method for DEG identification differed that in the previous study, which used a case-control method [15,16]. DEGs in the 0-8 h were mainly associated with virus adhesion, DEGs in the 8-24 h were largely associated virus endocytosis and transmission, and DEGs in the 24-72 h were involved in cell death and apoptosis. Therefore, the DEGs identified in the present study appear to be representative and accurate.
The CPE appeared as early as 8 h after CIK cells infected with GCRV. DEGs from this stage were mainly enriched with Focal adhesion and ECM-receptor interaction terms. Focal adhesions are sub-cellular structures that mediate the regulatory effects of cells in response to ECM adhesion [17]. In fact, previous studies showed that focal adhesion kinase can interact with rabies virus phosphoprotein P, which is involved in the process of viral infection [18]. During influenza A invasion of cells, inhibition of focal adhesion kinase signaling leads to reduced viral replication [19], suggesting that focal adhesion may be an important pathway during GCRV infection. Moreover, hspg2, ptpε, and ppl were significantly up-regulated during the 0-8 h post infection. HSPG2 plays an important role in mediating virus envelope-target cell interaction during human papillomavirus (HPV) and hepatitis C virus (HCV) infection [20,21]. PTPε plays an important role in the control of macrophage function [22]. And what's more, tyrosine-protein phosphatase non-receptor type 22 (PTPN22) associates with TNF receptor-associated factor 3 (TRAF3) to augment Toll-like receptor (TLR)-induced type I IFN production [23]. However, PPL acts as a cytolinker between intermediate filament scaffolding and the desmosomal plaque and is involved in epithelial cohesion, intracellular signal transduction and antigen presentation [24]. These results suggest that hspg2 may be involved in cellular binding of GCRV, after which GCRV likely binds to the receptor on the cell membrane, ptpε and ppl maybe as signaling molecules that regulate cellular processes.
During the 8 to 24 hpi period, enriched DEGs were mainly related to phagosome and lysosome pathways. In recent years, more and more studies have shown that the viral infection process is closely related to autophagy of host cells. Autophagy can be used as innate immunity and an adaptive immune response against intracellular pathogens [25]. In the phagosome pathway, many genes are up-regulated, including itgb2, srb1, and sec22b, and it is worth noting that the itgb2 is up-regulated at this period. Cell integrins are commonly used receptors for diverse viral pathogens [26], and reovirus internalization is mediated by integrin β1 proteins, most likely via clathrin-dependent endocytosis [27]. Moreover, integrin β3 proteins serve as receptors of pathogenic hantaviruses [28]. In the phagosome pathway, other receptor genes appear to be up-regulated, such as scavenger receptor srb1, which, along with CD81, is the receptor used for entry of the hepatitis C virus into liver cells [29]. The v-SNARE protein SEC22 is involved in the transport of eukaryotic cell endoplasmic reticulum (ER) and Golgi apparatus [30]. Hepatitis C virus may be assembled in the ER and transported to the Golgi apparatus in COPII vesicles to begin Golgi secretion [31]. However, a recent study shows that SEC22 reduces RNAi efficiency by influencing late endosomal functions, through the promotion of fusion between late endosomes and lysosomes [32]. These results suggest that itgb2 and/or srb1 may be important for cell entry of GCRV. Furthermore, GCRV entry into the endosomal region followed by fusion with lysosomes is likely to be regulated by sec22b.
In the 24-72 h post infection, DEGs were mainly enriched in lipid metabolism, programmed cell death, apoptosis, and inflammation. Genes in the cell death pathway have been studied intensely, and it is interesting to note that genes in the lipid metabolism pathway, such as cyp51a, are enriched. CYP51A1 is a cytochrome P450 enzyme, a conserved group of proteins that serve as key players in the metabolism of organic substances and the biosynthesis of important steroids, lipids, and vitamins in eukaryotes [33], and it plays an essential role in mediating membrane permeability [34]. These sterols, located on the plasma membrane of cells, play an important structural role in regulating the fluidity and permeability of membranes while also affecting the activity of enzymes, ion channels and other cellular components [35]. Therefore, these results suggest that GCRV probably not only causes cell death through cell apoptosis and inflammatory pathways but also changes the permeability of the cell membrane by affecting the structure of steroids in the cell membrane, leading to cell death through this mechanism.

CCK-8 Assay
CCK-8 detection kit (Beyotime, Shanghai, China) was used to measure cell viability according to the manufacturer's instructions. About 4 × 10 3 CIK cells were seeded in 96 well plates and cultured in M199 supplemented with 10% FBS at 28 • C for 24 h. Wells containing no cells but only culture medium alone were served as blank control. CIK cells were mock-infected or infected with GCRV at an MOI of 1 and maintained for 0, 8, 24, and 72 h at 28 • C, respectively. Then each well was added 10 µL of CCK-8 solution. After 4 h incubation at 28 • C, the plates were measured spectrophotometrically on Microplate Manager 6 (BIO-RAD, Hercules, CA, USA) at wavelength 450 nm. The relative cell viability was calculated as Ainfected-Ablank-i/A mock-Ablank % (Ablank-i represented the wells containing no cells but only medium and GCRV). Data were presented as means of at least three independent experiments ± standard deviation (S.D).

RNA Isolation, Library Construction, and Sequencing
Cells were collected at 0, 8, 24, and 72 h after infection (C0, C8, C24, and C72 groups, respectively), rapidly frozen in liquid nitrogen, and total RNA was extracted with Trizol reagent (Invitrogen, Waltham, MA, USA). RNase-free DNase treatment was performed in order to remove possible genomic DNA. The DNase was removed by phenol-chloroform re-extraction. mRNA was purified from total RNA using poly-T oligo-attached magnetic beads and segregated into 200-300 bp fragments. The concentration of mRNA used to prepare first strand cDNA is approximate 40 ng/µL and a total of 3 µg mRNA used as input material for first strand DNA synthesis. First strand DNA was synthesized from 6-base random primers and reverse transcriptase using RNA as a template, and second strand cDNA was synthesized using the first strand cDNA as template. In the second strand cDNA synthesis, the base T was replaced by U to generate a chain-specific library. RNA-Seq was performed using chain-specific kits, and chain-specific libraries were used to determine the transcriptional direction of sense and antisense strands to increase the accuracy of subsequent gene function annotation and gene expression analysis.
After library construction, library fragments were enriched by PCR amplification, and selection based on fragment size resulted in a 300-400 bp library. The library was tested using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA, USA), and both the total and effective concentration of the library were tested. Based on the effective concentration and the amount of data required, libraries containing different Index sequences (samples plus different Index sequences, and difference data of samples according to Index) were mixed proportionally. The mixed library was diluted to 2 nM and alkaline degeneration was performed to generate a single-chain library. Next, 150 bp pair-end reads were performed on the library using an Illumina NextSeq500 sequencing platform to carry out Next-Generation Sequencing (NGS) at Personalbio (Personalbio, Shanghai, China). For each treatment in each time point, the sample was sequenced three times as three technical replicates.

Data Analysis
Clean data were obtained by removing adapters, poly-N sequences and poor-quality data [36]. The Q20, Q30, and GC content of clean data were calculated, and all downstream analysis was performed using high-quality clean data.
Clean data were mapped to the grass carp reference genome [37] using TopHat2 software [38]. The number of reads mapped to each gene was used HTSeq software [39], and the reads per kilobase of the exon model per million mapped reads (RPKM) was calculated as expression values for each gene [40].

Differential Expression Analysis
Differential expression analysis using data collected at different time points was performed using the DESeq package [41]. The resulting p-values were adjusted using the Benjamini and Hochberg approach to control the false discovery rate. Genes with an adjusted p-value < 0.05 (q-value < 0.05) and fold change > 2 in DESeq analysis were assigned as differentially expressed genes (DEGs).
ClueGO and CluePedia were used to perform Gene Ontology (GO) enrichment analysis [42,43]. In GO enrichment analysis, only categories with a p-value < 0.05 were considered as enriched in the network.
The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to indicate the location of the DEGs in the different pathways. In this study, KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways [44]. Pathways with p-values < 0.05, were considered for as statistical significance.

Validation of DEGs by RT-qPCR
To confirm the reliability of the data obtained with RNA-Seq, eight DEGs were selected for RT-qPCR using the primers listed in Table S5. Primers were designed using the Primer Premier 5.0 software (PREMIER Biosoft International, Palo Alto, CA, USA) and only with an efficiency of 90-110% were used. First strand cDNAs were obtained from total RNA (>3 µg) using a random hexamer primer and the ReverTra Ace kit (TOYOBO, Osaka, Japan). RT-qPCR was performed using a fluorescence quantitative PCR instrument (BIO-RAD). Each RT-qPCR mixture contained 0.8 µL of each primer, 1 µL template, 10 µL 2× SYBR Green master mix (TOYOBO), and 7.4 µL ddH 2 O. The program for RT-qPCR was 95 • C for 10 s, followed by 40 cycles at 95 • C for 15 s, 60 • C for 15 s, and 72 • C for 30 s. Three replicates were included for each sample, and the β-actin gene of grass carp was served as an internal control to normalize the expression levels. Relative expression levels of the different genes were calculated using the 2 −∆∆Ct method [45]. All data represent the mean ± standard deviation of three replicates.

Conclusions
In conclusion, the results of this study provide insight into the mechanisms of grass carp reovirus infection. In the early stage of infection, GCRV interacts with polysaccharides and accumulates on the cell surface, then binds to the GCRV receptor to initiate endocytosis mediated by integrins. In the middle stage of infection, GCRV is likely to be regulated by v-SNARE protein and transported from the endosomal region to lysosomes. In the last stage, GCRV functions through several different pathways, by affecting the structure of steroids in the cell membrane, then leading to cell death in the final stage. These findings broaden our understanding of the invasion of GCRV and may assist the development of strategies to combat diseases caused by GCRV.