key: cord-0002934-n952re94
authors: Niu, Xiaosai; Wang, Yuyang; Li, Min; Zhang, Xiaorong; Wu, Yantao
title: Transcriptome analysis of avian reovirus-mediated changes in gene expression of normal chicken fibroblast DF-1 cells
date: 2017-11-25
journal: BMC Genomics
DOI: 10.1186/s12864-017-4310-5
sha: 602fbb4e85804b8ca5c0842b46afd343981181e4
doc_id: 2934
cord_uid: n952re94

BACKGROUND: Avian reovirus (ARV) is an important poultry pathogen that can cause immunosuppression. In this study, RNA-Seq technology was applied to investigate the transcriptome-wide changes of DF-1 cells upon ARV infection at the middle stage. RESULTS: Total RNA of ARV-infected or mock-infected samples at 10 and 18 h post infection (hpi) was extracted to build RNA-Seq datasets. Analysis of the sequencing data revealed that the expressions of numerous genes were altered, and a panel of differentially expressed genes were confirmed with RT-qPCR. At 10 hpi, 104 genes were down-regulated and 64 were up-regulated, while the expressions of 47 genes were increased and only one was down-regulated, which may play a role in retinoic acid biosynthesis, at 18 hpi in the ARV-infected cells. The similar profiles of up-regulated genes between the two groups of infected cells suggest that ARV infection activated a prolonged antiviral response of host cells. Alternative splicing analysis found no significantly changed events altered by ARV infection. CONCLUSIONS: Overall, the differential expression profile presented in this study can be used to expand our understanding of the comprehensive interactions between ARV and the host cells, and may be helpful for us to reveal the pathogenic mechanism on the molecular level. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (10.1186/s12864-017-4310-5) contains supplementary material, which is available to authorized users.

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 [1] . 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] [5] [6] . Efficient and simple detecting methods may be helpful to control and prevent ARV infection [7] . Horizontal transmission is the main route of infection, with infrequent egg transmission [8] . Though ARV was found to be ubiquitous in poultry flocks, several strains could cause severe diseases [9] . These pathogenic strains can cause tenosynovitis individually [9] , 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 [8] . 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 [15] . Another study demonstrated that ARV can induce autophagy to promote viral titer [16] . 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] [20] [21] [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] [28] [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.

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% CO 2 . 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 [16] . 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 (TCID 50 ) per milliliter (ml) of the virus was determined by the Reed-Muench method in CEF cells [30] .

One day before virus inoculation, approximately 2 × 10 6 DF-1 cells were seeded into 75 cm 2 flasks (Corning, ME, USA). When monolayer was complete (approximately 7 × 10 6 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 mockinfected 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) .

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.

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.

The index of the chicken reference genome (Ensembl, Galgal4, updated 11-2015) was built using Bowtie v2.2.3 [31] and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12 [32] . And HTSeq v0.6.1 was used to count the number of reads mapped to each gene [33] . Then, the expression level of each gene was calculated by the expected Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) [34] .

Differential expression analysis of two groups was performed using the DESeq R package (1.18.0) [35] . 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 [36] . 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 [37] .

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was carried out based on the basic rules of the MIQE guidelines [38] . 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 [39] .

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 [40] . Alternative splicing (AS) events were classified to five major types by the software rMATS (Multivariate Analysis of Transcript Splicing) v3.2.5 [41] . The number of AS events in each sample was estimated separately. Because the chicken genome has been recently updated [42] , the differentially expressed novel transcripts were retrieved by the BLAST tool on National Center for Biotechnology Information (NCBI).

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 Fig. 1 Overview of RNA-Seq approach. a Experimental setup for RNA-Seq datasets. b DF-1 cells were infected with ARV and the cytopathic effect was assessed at different time points. c The replication of ARV was monitored by RT-qPCR analysis of ARV genes M3 and S1 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 [45] . 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 (R 2 ), 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 [8] , there are no reads can be used to indicate the replication of ARV.

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 Table 2 with UniProtKB Keywords annotation [46] , as a 2-fold threshold is commonly used to indicate biological significance [47] . 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 [7] . 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 [43] . 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 [49] . 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 domaincontaining 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 acidinducible 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] [53] [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 IFNinduced 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 [55] . 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 [56] . Interestingly, another eukaryotic translation initiation factor 2-alpha kinase, EIF2AK3, was found to be downregulated 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 [57] . 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 [8] . 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 [60] . Later, studies demonstrated that the Wnt/β-catenin signaling pathway is necessary and sufficient to induce early steps of synovial joint formation [61] . Subsequently, a precise expression pattern of various Wnts was analyzed during chick wing development [62] . 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 [60] . ARV can replicate and perhaps be persistent at hock joint of chicken [48] . The up-regulation of Wnt14, combined with the induction of apoptosis [15] , may be responsible for ARVinduced joint damage and more severe tendon rupture.

The raw data of the RNA-Seq datasets was deposited in Sequence Read Archive (SRA) under the accession number SRP104704. 

Virus taxonomy: ninth report of the international committee on taxonomy of viruses

The history of avian reovirus

Viral arthritis/tenosynovitis: a review

Occurrence of reovirus infection in Muscovy ducks (Cairina Moschata) in south western Poland

Phenotypic, genotypic and antigenic characterization of emerging avian reoviruses isolated from clinical cases of arthritis in broilers in Saskatchewan

Flock-level prevalence, geographical distribution, and seasonal variation of avian reovirus among broiler flocks in Ontario

The development and evaluation of cross-priming amplification (CPA) for the detection of avian reovirus (ARV)

Avian immunosuppressive diseases and immunoevasion

Avian reovirus: structure and biology

Blue wing disease of chickens: isolation of avian reovirus and chicken anaemia agent

Blue wing disease of chickens: experimental infection with a Swedish isolate of chicken anaemia agent and an avian reovirus

Replication of four antigenic types of avian reovirus in subpopulations of chicken leukocytes

Avian reovirus σC protein induces apoptosis in cultured cells

Extensive syncytium formation mediated by the reovirus FAST proteins triggers apoptosis-induced membrane instability

Avian reovirusinduced apoptosis related to tissue injury

Reovirus triggers autophagy in primary chicken fibroblast cells and Vero cells to promote virus production

RhoA/ROCK1 regulates avian Reovirus S1133-induced switch from autophagy to apoptosis

Autophagy inhibitors reduce avian-reovirus-mediated apoptosis in cultured cells and in chicken embryos

Avian reovirus influences phosphorylation of several factors involved in host protein translation including eukaryotic translation elongation factor 2 (eEF2) in Vero cells

Retardation of cell growth by avian reovirus p17 through the activation of p53 pathway

Cell entry of avian reovirus follows a caveolin-1-mediated and dynamin-2-dependent endocytic pathway that requires activation of p38 mitogen-activated protein kinase (MAPK) and Src signaling pathways as well as microtubules and small GTPase Rab5 protein

The p17 nonstructural protein of avian reovirus triggers autophagy enhancing virus replication via activation of phosphatase and tensin deleted on chromosome 10 (PTEN) and AMP-activated protein kinase (AMPK), as well as dsRNA-dependent protein kinase (PKR)

Avian reovirus S1133-induced DNA damage signaling and subsequent apoptosis in cultured cells and in chickens

Proteomics analysis of the DF-1 chicken fibroblasts infected with avian reovirus strain S1133

RNA-seq and microarray complement each other in transcriptome profiling

Transcriptome analysis of Pinus Monticola primary needles by RNA-seq provides novel insight into host resistance to Cronartium Ribicola

RNA-Seq analysis reveals new gene models and alternative splicing in the fungal pathogen Fusarium graminearum

TRPM2 SNP genotype previously associated with susceptibility to Rhodococcus equi pneumonia in quarter horse foals displays differential gene expression identified using RNA-Seq

Using galaxy-P to leverage RNA-Seq for the discovery of novel protein variations

A simple method of estimating fifty percent endpoints

Fast gapped-read alignment with bowtie 2

TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions

HTSeq-a python framework to work with high-throughput sequencing data

Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms

Differential expression analysis for sequence count data

Controlling the false discovery rate: a practical and powerful approach to multiple testing

2.0: a web server for annotation and identification of enriched pathways and diseases

The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments

Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method

Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks

rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data

Differential expression profile of chicken embryo fibroblast DF-1 cells infected with cell-adapted infectious Bursal disease virus

Host cell transcriptome profile during wild-type and attenuated dengue virus infection

Immune-Related Gene Expression Patterns in GPV-or H9N2-Infected Goose Spleens

the SIB Swiss Institute of Bioinformatics and the Protein Information Resource (PIR)

Transcriptome-wide analysis of hepatitis B virus-mediated changes to normal Hepatocyte gene expression

Avian reovirus infection

EPSTI1 is involved in IL-28A-mediated inhibition of HCV infection

Activation of the TLR3 pathway regulates IFNbeta production in chickens

Polymorphisms of chicken TLR3 and 7 in different breeds

MDA5 can be exploited as efficacious genetic adjuvant for DNA vaccination against lethal H5N1 influenza virus infection in chickens

Chicken cells sense influenza a virus infection through MDA5 and CARDIF signaling involving LGP2

Mediates activation of the IFN gene independently of the RIG-I gene

ERK2 and RSK2 mediate phosphorylation of PKR (Thr451) and PKR directly catalyzes the phosphorylation of eIF2α at Ser51

Possible involvement of the double-stranded RNA-binding Core protein ςA in the resistance of avian Reovirus to interferon

Regulated translation initiation controls stress-induced gene expression in mammalian cells

Protein synthesis and endoplasmic reticulum stress can be modulated by the hepatitis C virus envelope protein E2 through the eukaryotic initiation factor 2α Kinase PERK

The SARS Coronavirus 3a protein causes endoplasmic reticulum stress and induces ligand-independent downregulation of the type 1 interferon receptor

Wnt-14 plays a pivotal role in inducing synovial joint formation in the developing appendicular skeleton

Wnt/beta-catenin signaling is sufficient and necessary for synovial joint formation

Comparative analysis of the expression patterns of Wnts during chick limb development

Thanks to Huaisheng Wu for the assistance in performing RNA-Seq analysis. 

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. Ethics approval and consent to participate Not applicable.

Not applicable.

The authors declare that they have no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Submit your next manuscript to BioMed Central and we will help you at every step: