Avian reovirus (ARV) is member of the Orthoreovirus genus that has recently been classified into the Spinareovirinae subfamily, which is one of two subfamilies in the Reoviridae family. ARV is an important pathogen of birds and has been impacting poultry for nearly 60 years since it was first detected in 1957 [2, 3], and it is still prevalent in poultry until now, causing considerable economical loss in the global poultry industry [4–6]. Efficient and simple detecting methods may be helpful to control and prevent ARV infection. Horizontal transmission is the main route of infection, with infrequent egg transmission. Though ARV was found to be ubiquitous in poultry flocks, several strains could cause severe diseases. These pathogenic strains can cause tenosynovitis individually, and additionally usually cause mixed infections together with other pathogens, such as chicken anemia agent [10, 11]. It has been demonstrated that ARV can replicate in macrophages and cause immunosuppression [8, 12].
The pathogenicity and epizootiology of ARV have been well studied, but the pathogenesis at the molecular level is poorly understood. An excellent review on the structural and biological characteristics of ARV was published 10 years ago. Though many researchers have done brilliant work to reveal the pathogenesis of ARV infection at the molecular level in recent years, several major questions raised in the review remain unresolved. Previous studies showed that σC and P10 can induce apoptosis in different ways [13, 14], and a subsequent study correlated ARV-induced apoptosis with tissue injury. Another study demonstrated that ARV can induce autophagy to promote viral titer. Subsequent studies revealed the connection between ARV-induced autophagy and apoptosis [17, 18]. It was also demonstrated that ARV disrupts many cellular pathways, regulating protein translation, cell proliferation, and cell metabolism [19–22]. However, these results are scattered and hard to reconcile. Some studies applied proteomic analysis and microarray analysis to get a comparatively integrated data set [23, 24]. However, these methods have several disadvantages compared with RNA-Seq. RNA-Seq now provides a way to investigate virus-mediated changes on the transcriptome of host cells, with high accuracy and low background [25, 26]. Additionally, it provides information on alternative splicing events, analyzing single nucleotide polymorphism, and predicting novel transcripts [27–29]. In this study, we tried to build a complete expression profile of ARV-mediated changes at the transcriptional level using RNA-Seq to unveil the complex interactions between ARV and host cells.
Cell culture and virus inoculation
Chicken embryonic fibroblast cell line DF-1 (CRL-12203, ATCC) cells were cultured with high glucose (4.5 g D-Glucose/L) Dulbecco’s Modified Eagle Medium (HG-DMEM) (Basal Media, Shanghai, China) supplemented with 10% (v/v) fetal bovine serum (FBS) (Gibco, Shanghai, China) at 37 °C and 5% CO2. The ARV strain GX/2010/1, causing severe tenosynovitis and enteritis, was isolated by our lab and propagated in chicken embryo fibroblasts (CEF) cells, and was reported to trigger autophagy in host cells to promote virus production. The sequences of this strain are available in the GenBank database under the accession numbers KJ476699−KJ476708. The sixth generation of the purified virus was used in this study and the median tissue culture infective dose (TCID50) per milliliter (ml) of the virus was determined by the Reed-Muench method in CEF cells.
One day before virus inoculation, approximately 2 × 106 DF-1 cells were seeded into 75 cm2 flasks (Corning, ME, USA). When monolayer was complete (approximately 7 × 106 cells), culture medium was discarded and the cells were rinsed with phosphate buffered saline once. The purified virus was diluted to 10 multiplicity of infection (MOI) per 5 mL with HG-DMEM and applied into each flask of the ARV-infected group and an equal volume of HG-DMEM was added in the mock-infected group. After incubated at 37 °C for 1.5 h, the medium was changed to HG-DMEM supplemented with 2% FBS. Then the cells were continued to be incubated at 37 °C for 2, 10, 18 and 24 h (Fig. 1a).
Total RNA extraction and cDNA library construction
At specified hours post infection (hpi), the medium was discarded and total RNA was extracted from triplicate samples of uninfected or ARV infected groups using the Ultrapure RNA Kit (CWBIO, Beijing, China) according to manufacturer’s protocol. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample.
Clustering and sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 150 bp paired-end reads were generated. Raw reads of fastq format were first processed through custom written Perl scripts. At the same time, Q20, Q30, and GC content were calculated. All of the downstream analyses were based on clean, high quality data.
Reads mapping and quantification of gene expression level
The index of the chicken reference genome (Ensembl, Galgal4, updated 11-2015) was built using Bowtie v2.2.3 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12. And HTSeq v0.6.1 was used to count the number of reads mapped to each gene. Then, the expression level of each gene was calculated by the expected Fragments Per Kilobase of transcript per Million fragments mapped (FPKM).
Differential expression analysis
Differential expression analysis of two groups was performed using the DESeq R package (1.18.0). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on a negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq were considered to be differentially expressed. Additionally, KOBAS 2.0 software was used to test the statistical enrichment of differentially expressed genes (DEGs) in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.
Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was carried out based on the basic rules of the MIQE guidelines. Briefly, 5 μg of total RNA (described above) was reverse transcribed using M-MLV reverse transcriptase (Transgen, Beijing, China) with a random hexamer primer (Genscript, Nanjing, China). The mixtures were diluted 1:10 with nuclease free water and then used as templates for qPCR. The qPCR analysis was performed using AceQ® qPCR SYBR® Green Master Mix (Vazyme, Nanjing, China) with 250 nM forward and reverse primers (Additional file 1). The reaction was carried out using LightCycler® Nano (Roche) with the following cycling conditions: an initial denaturation at 95 °C for 600 s followed by 45 cycles of 95 °C for 10 s and 60 °C for 30 s. Fold change was determined by the 2-△△Ct method.
Novel transcripts and alternative splicing prediction
The Cufflinks v2.1.1 Reference Annotation Based Transcript (RABT) assembly method was used to construct and identify both known and novel transcripts from TopHat alignment results. Alternative splicing (AS) events were classified to five major types by the software rMATS (Multivariate Analysis of Transcript Splicing) v3.2.5. The number of AS events in each sample was estimated separately. Because the chicken genome has been recently updated, the differentially expressed novel transcripts were retrieved by the BLAST tool on National Center for Biotechnology Information (NCBI).
ARV infection of DF-1 cells and viral replication dynamics
To further study the molecular mechanism of ARV infection, DF-1 cells were infected with the virus for different time points at 10 MOI. The high dosage of virus was used to overcome the influence of uninfected cells [43, 44]. Infection and mock-infection were performed in biological triplicate for each time point and total RNA was extracted from both groups. The replication of the viral genome was determined by RT-qPCR and the fold change of M3 and S1 showed similar trends (Fig. 1c). Cytopathic effects could be seen at 18 hpi (Fig. 1b). To obtain an obviously changed transcriptome profile and minimize the influence of cell death and lysis, data at 2 hpi and 24 hpi were discarded and the remaining two groups were analyzed by RNA-Seq. One uninfected 10 hpi sample was lost, leaving a final total of 11 samples that were sequenced.
After an overall quality review, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads, and then the cDNA library was constructed, with quality assessment. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 150 bp paired-end reads were generated. The sequence run of each sample yielded at least 42 million clean reads and the lowest value of the reads possessing a Q-score > 20 was 95%, and the bottom line of the reads with a Q-score > 30 was 88% (Table 1). These results meet the requirements that more than 10 million reads are needed to construct a high quality eukaryotic transcriptome for discovering new genes and transcripts. Importantly, all samples had between 71.81% and 77.45% of total reads mapped to the chicken reference genome, and the percentage of uniquely mapped reads was between 70.75% and 76.26% (Table 1).
Then, the mapped data were used to predict the novel transcripts and analyze the five major types of AS events, including SE (skipped exon), A5SS (alternative 5′ splice site), A3SS (alternative 3′ splice site), MXE (mutually exclusive exon), and RI (retained intron). The predicted novel genes were further analyzed together with the known genes. No significantly changed AS events were found between each group. The mapped data was normalized by calculating the FPKM and the distribution of mean FPKM per gene was found to be uniform between the four conditions (Fig. 2a). The correlation of gene expression levels between all of the samples was investigated using the squared Pearson correlation coefficient (R2), and the minimum value was 0.979 (Fig. 2b). These results indicate that the expression levels of different genes or groups of genes are comparable, suggesting that the treatment is repeatable and has little variation. Therefore, the accuracy of the subsequent analysis of differentially expressed genes is likely to be high. Because the transcripts of the ARV genome do not have a poly-A 3′ tail, there are no reads can be used to indicate the replication of ARV.
Differentially expressed genes upon ARV infection
To further investigate the differential expression patterns in DF-1 cells between infected and mock-infected samples, the normalized gene expression level data were analyzed by DESeq. The resulting P-values were adjusted after correction for multiple testing and the DEGs were defined by having adjusted P-values (padj) < 0.05. The infected and mock-infected were compared with each other and the outline of the DEGs can be seen in Fig. 3a. Though there were significant alterations in the pairwise comparisons of different time points in the infected samples (Fig. 3b), the corresponding mock-infected samples also had big differences (Fig. 3c). This interference might result from cell culture and so should be discarded in the future studies. All of the DEGs were clustered and the results exhibited a clear time-dependent change in gene expression (Additional file 2).
The distinct effect on gene expression upon ARV infection was carefully examined. Compared with mock infected controls, 168 changes in the transcriptome, with 64 up-regulated and 104 down-regulated DEGs, were observed in response to ARV infection at 10 hpi (Fig. 3d, Additional file 3). Interestingly, only 47 up-regulated DEGs, with 31 genes in accordance with the 10 hpi group, and a single down-regulated novel gene were identified at 18 hpi (Fig. 3e, Additional file 3). These novel genes among the DEGs were retrieved by the BLAST tool. In addition, the pathway enrichment result can be seen at additional file 4. The DEGs with a fold change larger than 2 are listed in Table 2 with UniProtKB Keywords annotation, as a 2-fold threshold is commonly used to indicate biological significance.
Four up-regulated DEGs and two down-regulated DEGs were selected to be validated with RT-qPCR. The results show a similar pattern of ARV-mediated changes as was seen in the DEG analysis of RNA-Seq data (Fig. 4).
ARV is one of the major pathogens that can cause immunosuppression in poultry. Though the pathology and some molecular characteristics of ARV have been well studied [2, 8, 48], there are only a few reports that can be used to help us understand the molecular basis of ARV infection. In this study, the DEG listed in Table 2 show that DF-1 cells exerted a prolonged antiviral response upon ARV infection. The up-regulation of RSAD2 (radical S-adenosyl methionine domain containing 2), IFIT5 (interferon induced protein with tetratricopeptide repeats 5), OASL (2′-5′-oligoadenylate synthetase-like), ISG12(2) (interferon-stimulated genes) and Mx (myxovirus resistance) have been reported in the infection of infectious bursal disease virus (IBDV), which is another important pathogen similar to ARV but can cause much higher mortality and much more serious immunosuppression. EPSTI1 (epithelial-stromal interaction 1) does not have a Gene Ontology (GO) annotation. A recent report indicated that EPSTI1 plays a key role in IL-28A (interferon-λ2) mediated antiviral activity. These interferon- (IFN-) induced genes (ISGs) with high expression levels reflect the stimulation of IFNs. Even though no significant elevation of expression levels of IFN genes were identified in this study, the up-regulation of TLR3 (Toll-like receptor 3), MYD88 (Myeloid differentiation factor 88), IRF1 (IFN regulatory factor 1), and IRF3 (IFN regulatory factor 3) were found. TLR3 plays key roles in detecting virus-derived dsRNA and the TLR3 genes are polymorphic among different chicken breeds [50, 51]. In addition to TLR-induced pathways, members of the RLR family (retinoic acid inducible gene-I like receptor) constitute another TLR-independent anti-virus system. In our results, DHX58 (DEXH-box helicase 58, also known as LGP2, laboratory of genetics and physiology 2, or RLR3, RIG-I like receptor 3), IFIH1 (IFN-induced helicase C domain-containing protein 1, also known as MDA5, melanoma differentiation-associated protein 5), TRIM25 (Tripartite motif-containing protein 25), and a predicted TRIM25-like gene were found to be up-regulated at 10 hpi or at both 10 and 18 hpi time points. Chickens lack RIG-I (retinoic acid-inducible gene I), but the function of sensing viral infections can be performed by LGP2 and MDA5, which can interact with MAVS (mitochondrial antiviral signaling protein) or STING (stimulator of IFN genes) to stimulate the expression of IFNs [52–54].
DF-1 cells construct an antiviral environment through the expression of ISGs, including EIF2AK2 (Eukaryotic translation initiation factor 2-alpha kinase 2, also known as PKR, protein kinase RNA-activated). PKR is IFN-induced dsRNA-dependent enzymes. Active PKR can catalyze Ser-51 phosphorylation of the alpha subunit of EIF2, resulting in inhibition of protein synthesis at the initiation step of translation. However, a previous report demonstrated that σA, an ARV encoded dsRNA binding protein, can block the activation of PKR and restore translation. In that report, the inhibition of vaccinia virus replication might reflect mechanisms other than OAS and PKR to be responsible for the antiviral effects. Interestingly, another eukaryotic translation initiation factor 2-alpha kinase, EIF2AK3, was found to be down-regulated in our results. EIF2AK3, also known as PERK (PKR-like endoplasmic reticulum-resident kinase), is one of eIF2α kinases regulating gene expression in the unfolded protein response (UPR) and in amino acid starved cells. Protein synthesis can be inhibited during viral infection due to ER stress triggered by UPR, and different viruses may adapt different strategies to interfere with the activity of PERK [58, 59]. The depressed expression of PERK may reflect that ARV can impair the stress response and activate protein translation in DF-1 cells. The regulation of the host cell translation system ensures efficient replication of ARV. There is also a hypothesis that gene expression of ARV is mainly regulated at the translational level, rather than transcriptional level. The replication level of viral genome determined by RT-qPCR in our results is consistent with this hypothesis.
ARV was initially detected from the clinical case of tenosynovitis, and a direct link between the virus and disease had been conclusively demonstrated [2, 8]. Though this virus has been studied for many years, the molecular pathogenesis of the disease remains unclear. In our results, the elevated expression of a gene, WNT9a (also known as Wnt14), was observed, which might play a key role in the development of the disease. A previous report identified that Wnt14 plays a pivotal role in initiating synovial joint formation in the chick limb, but the researchers were unable to determine the specific pathway that is responsible for transducing the Wnt14 signal in joint formation. Later, studies demonstrated that the Wnt/β-catenin signaling pathway is necessary and sufficient to induce early steps of synovial joint formation. Subsequently, a precise expression pattern of various Wnts was analyzed during chick wing development. Continued expression of Wnt14 in the mature joint might be good for the maintenance of joint integrity and was presumed to play a role in the etiology of rheumatoid arthritis in humans. ARV can replicate and perhaps be persistent at hock joint of chicken. The up-regulation of Wnt14, combined with the induction of apoptosis, may be responsible for ARV-induced joint damage and more severe tendon rupture.
In conclusion, our results show that ARV infection stimulates a prolonged antiviral response in host cells and interferes with cell growth and cell death pathways. Our results also provide information that may be helpful to further investigate the pathogenesis of ARV infection. Combined with previous studies, we can begin to piece together the interactions between ARV and host cells (Fig. 5). However, the details of these interactions need to be further investigated in future studies.