Differential Expression Profile of Chicken Embryo Fibroblast DF-1 Cells Infected with Cell-Adapted Infectious Bursal Disease Virus

RNA-Seq was used to unveil the transcriptional profile of DF-1 cells at the early stage of caIBDV infection. Total RNAs were extracted from virus-infected cells at 0, 6 and 12 hpi. RNA-Seq datasets of respective samples mapped to 56.5–57.6% of isoforms in the reference genome Galgal4.73. At 6 hpi, 23 isoforms underwent an elevated expression, while 128 isoforms were up-regulated and 5 were down-regulated at 12 hpi in the virus-infected group. Besides, 10 isoforms were exclusively expressed in the virus-infected cells. Though no significant change was detected in cytokine and interferon expression levels at the first 12 hours of infection, modulations of the upstream regulators were observed. In addition to the reported regulatory factors including EIF2AK2, MX, OAS*A, GBP7 and IFIT, IBDV infection also triggered a IFIT5-IRF1/3-RSAD5 pathway in the DF-1 cells which potentially restricted the viral replication cycle in the early infection stage. Over-expression of LIPA and CH25H, together with the suppression of STARD4, LSS and AACS genes implied a modulation of membrane fluidity and lipid raft arrangement in the infected cells. Alternative splicing of the EFR3 homolog A gene was also through to be involved in the lipid membrane regulation, and these cumulative responses projected an inhibition of viral endocytosis. Recognition of viral RNA genomes and intermediates was presumably enhanced by the elevated levels of IFIH1, DHX58 and TRIM25 genes which possess properties on detecting viral dsRNA. On the other hand, the caIBDV arrested the host's apoptotic process by inducing the expression of apoptosis inhibitors including NFKBIA/Z, TNFAIP2/3 and ITA at the first 12 hours of infection. In conclusion, the differential expression landscape demonstrated with RNA-Seq provides a comprehensive picture on the molecular interactions between host cells and virus at the early stage of infection.


Introduction
Infectious bursal disease (IBD) has been striking chicken flocks for more than fifty years exerting an considerable economical impact to the global poultry industry. The disease brings a direct mortality ratio up to 90-100% [1,2], and as it causes destruction of B-lymphocytes in the bursa of Fabricius, it leads into severe immunosuppression and hence secondary infections may result in infected chickens [3,4,5]. Infectious bursal disease virus (IBDV) is the causative agent of the disease. Two serotypes are identified in which serotype 1 comprises pathogenic strains, whereas serotype 2 strains cause neither disease nor protection against serotype 1 strains in chickens [6,7,8,9]. It is demonstrated that the virus propagates in the actively proliferating IgM-bearing B-lymphocytes and hence induces apoptotic effects [10,11,12]. Though the pathogenicity and epizootiology have been studied for a certain period of time, the molecular interactions between the host cells and the viruses have not been well defined yet. In recent years studies have started to focus on the molecular mechanisms involved in the host responses upon IBDV infection. Quantitative RT-PCR (qRT-PCR) and microarray assays are increasingly employed to reveal the transcriptional changes of the host cells in response to IBDV infections [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. While some studies also utilize proteomic approaches to identify the differentially expressed protein during the course of IBDV infection [31,32]. Majority of these studies emphasized the cytokine responses including interleukin and interferon expressions, whereas some of these studies revealed expression of mRNA related to apoptotic mechanisms. Up to now, however, there is no comprehensive transcriptional landscape described in the cells upon IBDV infection. In order to explore the differential expression pattern in the event of IBDV infection, RNA sequencing (RNA-Seq) was employed to assay the transcript variations across the entire chicken genome. RNA-Seq reveals a high overall sensitivity on differentially expressed gene level compared with other whole-transcriptome expression quantification platforms including microarrays [33,34]. The prerequisite of hybridization-based microarray assays relies on existing knowledge about genome sequences [35,36] and hence limits the detection of novel, rare transcript species exist in the transcriptome. Whereas RNA-Seq takes an advantage not only in determining the differential expression level of transcripts, but it also provides evidence on transcript splice-variants, isoforms and single nucleotide polymorphism (SNPs) [37]. It has also been demonstrated that RNA-Seq is highly accurate for determining gene expression levels as performed with qPCR [38]. Background levels resulting from cross-hybridization is also much lower than occurred in microarray assay [39]. Taking these advantages, in this study we made use of RNA-Seq to unveil the transcriptomic dynamics upon caIBDV infection in DF-1 cells and to reveal a more comprehensive molecular interactions between the host cells and the virus.

Virus inoculation
DF-1 cells were seeded into each well of 6-well plate (Costar 3516, Corning, NY) at 1 × 10 4 cells at 24 hrs prior to virus inoculation. A total of six individual wells were prepared for each sampling time point of both mock-and caIBDV-infected groups. The purified caIBDV was diluted to 10 multiplicity of infection (m.o.i.) per ml with DMEM-HG without serum supplement. Before inoculation, culture medium was discarded and the cells were rinsed with 1 × phospate-buffered saline once. One ml of diluted caIBDV was applied into each well of the virus-infected group, whereas DMEM-HG without serum supplement was added in the mockinfected group. At 0, 6, 12 hrs post-infection (hpi), medium was aspirated from the wells of mock-and caIBDV-infected groups. One ml of SV RNA Lysis Buffer (Promega, WI) was added into each well and the cells were dispersed and lysed with repeated pipetting. A total of six individual cell lysates were collected from each treatment group at designated time points and kept at -80°C for total RNA extraction thereafter.

Total RNA isolation
Total RNAs were extracted from cell lysates with SV96 Total RNA Isolation System (Promega, WI) according to manufacturer's protocol. Total RNAs extracted from each treatment group at the same time point were pooled into one tube for the downstream processes. A total of six RNA samples were obtained representing two treatments at three time points. Quantity of the pooled RNA samples was determined with Quant-iT RiboGreen RNA Reagent and Kit (Molecular Probes, OR), while the quality and integrity of RNA were confirmed with spectrophotometry and formaldehyde 1% agarose gel electrophoresis. All total RNA samples possessed A260/ A280 ratio > 1.8 and the 28S:18S rRNA intensity ratio of 1.76 estimated with Nanodrop 2000 (Thermo Scientific, DE) and Quantity One Analysis Software (Bio-Rad, CA) respectively. Whole transcriptome shotgun sequencing A 2 × 250 paired-end sequencing run was conducted using MiSeq Reagent Kit v2, 500 cycles with a standard flow cell (14 tile) (Illumina, CA) on an Illumina MiSeq machine. Raw reads generated from the run were demultiplexed and underwent quality trimming prior to the downstream analytical pipelines. Only reads with Q Score ! 30, and aligned with corresponding paired reads were retained for the differential expression analysis.

Sequence analysis
The paired-end sequence reads of individual sample were aligned and mapped to a chicken reference genome (Gallus Gallus reference genome Galgal4.73) with STAR RNA-Seq aligner 2.3.0 [41] in default setups. Both the unmasked reference genome and the annotations of both coding and non-coding genes were obtained from Ensembl database [42]. The mapped reads were then converted to binary BAM files with SAMtools 1.0.2 [43] for transcript expression levels comparison. The aligned, mapped reads from the same time point of both mock-and caIBDVinfected group were compared and the expression levels at the isoform and gene level were calculated with Cuffdiff package in Cufflink transcript assembly tool [44]. The same annotation and reference genome (Galgal4.73) files as for STAR aligner were used in Cuffdiff pipeline. Expression level of each transcript in the samples were expressed in fragments per kilobase of transcript per million mapped fragments, FPKM. Each observed alternation in transcript expression with p<0.05 was regarded as statistical significant. Quality of the libraries and sequencing performance were examined with RNASeQC tool [45]. The differentially expressed genes were grouped with hierarchical clustering algorithms [46] in the RNA-Seq analysis package of Gene Pattern web-based tool (Broad Institute, MIT).
The potential protein-protein inactions among the genes in the altered expression profiles derived from RNA-Seq analyses were predicted with STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) v9.05 [47].
In addition to the chicken reference genome, the paired-end sequence reads of each sample group were also mapped to a IBDV reference sequence (IBDV Segment A, NC_004178.1; IBDV Segment B, NC_004179.1 [48]) with GS Reference Mapper in 454 Sequencing System Software Package v2.7 (Roche, CT) in order to assess the change of caIBDV viral load in the infected DF-1 cells.

Nucleotide sequence accession number
The RNA-Seq datasets have been deposited at NCBI Gene Expression Omnibus (GEO) repository under the accession number of GSE60268.

Infection experiment and RNA-Seq
To study the transcriptomic profile of chicken embryonic fibroblast DF-1 cells during caIBDV infection, caIBDV strain D78 at 10 m.o.i. was used to infect DF-1 cells for 0, 6 and 12 hours. The relatively high dosage of virus was used for inoculation as to minimize the influence exerted from a large number of uninfected cells [49].
Total RNAs were isolated from both mock-and caIBDV-infected groups at respective time points and cDNA libraries were constructed for a 2 × 250 paired-end sequencing run on an Illumina MiSeq machine. The sequence run yielded a total of about 32 millions passed-filter reads, and 85.1% of the reads possessed Q-Score ! 30. In average, > 4 million quality-trimmed reads were generated for each sample (Table 1). Though it has been suggested that high-quality eukaryotic transcriptome reconstruction requires more than ten million reads for discovering new genes and transcripts [34], the relatively shallow RNA-Seq datasets are otherwise sufficient to provide accurate differential expression trends of transcript as confirmed with qRT-PCR [50,51,52]. Hence RNA-Seq datasets obtained in this study can illustrate the evidence on differential expression landscape between the mock-and virus-infected conditions.

Viral replication dynamics
After demulitplexing and quality trimming, the six datasets were mapped to IBDV reference genome (Segment A, NC_004178.1; IBDV Segment B, NC_004179.1) with GS Reference Mapper Algorithms (Roche, CT). The mapping matrix indicates that no sequence read from the mock-infected group matched to the reference genome. On the other hand, unique sequence reads of the segment A of the IBDV genome was detected in 6hpi and 12hpi, whereas those matched to segment B was observed at 12hpi (Fig 1). The sequence reads covered 20.01% to 63.87% of the entire segment A sequence at 6hpi and 12hpi samples respectively, while 56.17% of segment B sequence was retrieved from the 12hpi sample. This finding reveals that the replication machinery of caIBDV started as early as at 6 hours after the inoculation. The segment A sequences identified in the dataset did not only indicate the replication of IBDV genomic RNAs, but it also reflected the transcription of sub-genomic RNA species encoding for the viral structural proteins VP2, 3, 4 and 5 for the assembly of a mature virion. While the segment B sequences were observed at 12hpi, indicating that there was an expression of viral RNA polymerase VP1 and also the replication of this segment. The increase of mapped sequence reads indicates an active viral replication process in this 12-hours post inoculation period.

Differential expression patterns
To identify the differentially expressed transcripts in DF-1 cells upon caIBDV infection, the six RNA-Seq datasets were then aligned and mapped to a chicken reference genome Ensembl Gal-gal4.73 with STAR-aligner. The numbers of genes and isoforms mapped to the reference genome respective to each dataset remained constant throughout the course of experiment ( Table 2). The datasets mapped to 57.7-58.9% of genes and 56.5-57.6% of isoforms in the reference genome. There was no significant alternations in the expressed gene/ transcript isoform Besides, 9 additional isoforms were detected in the caIBDV-infected cells only (Fig 3). The result also demonstrates that there was a switch of alternatively spliced transcripts of EFR3A gene between mock-and caIBDV-infected groups (Fig 4).
The differentially expressed isoforms were classified into 12 clusters according to their expression patterns in hierarchical clustering (Fig 5). Early expressed transcripts (i.e., at 6 hpi) were observed in mainly in cluster 1, 2, 3 and 9, while cluster 12 comprised 5 down-regulated isoforms. Despite most of the expression patterns of these clusters are similar to each other and show limited co-relation with the potential biological process activated by the caIBDV-infection event according to their GO term accessions (Table 3), the result reveals that transcripts involved in anti-viral responses, transcription regulation and membrane proteins were dominantly expressed in caIBDV-infected cells at the very early stage of infection. While transcripts related to apoptotic activities, interferon-dependent immune responses and inflammatory responses were tuned up at 12 hpi.

Anti-viral responses exerted by DF-1 cells
It has been well demonstrated that host cells execute a cascade of machinery to combat against viral particle invasion. Previous studies have shown that IBDV infection triggered expression of cytokine genes that are typically up-regulated in general viral infections. Expression of T-      expressions including IL-4, IL-5, IL-13 was also reported [24]. Genes of cytokine that initiate inflammatory responses including IL-8, nitric oxide synthase (iNOS), and cyclooxygnease (COX-2) were also shown to be up-regulated [17,18,27]. In this study, however, no significant elevation of expression levels of these cytokine genes was identified in caIBDV-infected cells.
In contrary, the regulatory factors modulating these cytokines including the well characterized eukaryotic translation initiation factor 2-alpha kinase 2 (EIF2AK2, or protein kinase R PKR), interferon-induced GTP-binding protein Mx1 (MX), 59 kDa 2'-5'-oligoadenylate synthetaselike protein (OAS Ã A), guanylate binding protein (GBP7) and interferon-induced protein with tetratricopeptide (IFIT) were profoundly expressed. Among the 139 differentially expressed transcripts, IFIT5 (interferon-induced protein with tetratricopeptide repeats 5) was expressed in the highest fold change during caIBDV infection. IFIT5 is a member of IFIT1 family which expression can be induced by virus infection, interferons, dsRNA and lipopolysaccharides [53,54]. It was demonstrated that IFIT5 potentiates anti-viral responses by promoting the interferon regulatory factor 3 (IRF3)-or nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB)-mediated gene expressions [55]. Over-expression of IFIT5 in the infected DF-1 cells was observed accompany with the elevation of IRF3 at 6 hpi, and progressively with other IRF-3-mediated genes at 12 hpi (Fig 6). One of the dominant anti-viral components potentially triggered by ITIF5-IRF3 pathway was RSAD2 (radical SAM domain-containing 2, viperin) gene. RSAD2 gene was found expressed only in caIBDV-infected cells but not in mock-infected cells at 12 hpi (Fig 3). It is an endoplasmic reticulum-associated virus inhibitory protein which can be induced by type I, II and III IFNs, double-stranded (ds) DNA, dsRNA analog poly I:C, lipopolysaccharide (LPS) and infection of various viruses [56][57][58][59][60][61]. Besides, it was shown that RSAD2 could be induced via an IFN-independent pathway by transcription factor IRF-1 [62]. Definite anti-viral machineries exerted by RSAD2 has not been concluded yet, but it is believed that it involves in inhibiting viral replication indirectly by altering the cell survival control [63], and membrane fluidity modulation which prevents the budding of viruses [64]. So the switching-on of RSAD2, presumably by ITIF5-IRF3 or by IRF1 regulation, may contribute to the early viral defense in IBDV infection. Apart from the defense effectors, genes possessed with accessory anti-viral functions were also regulated in caIBDV-infected cells. The interferon-induced helicase C domain-containing protein 1 (IFIH1, or melanoma differentiation-associated protein 5, MDA5) was found expressed in 3.01 and 8.27 fold at 6 and 12 hpi respectively. It was reported that IFIH1 interacts with probable ATP-dependent RNA helicase (DHX58, or Laboratory of Genetics and Physiology 2, LGP2; expressed to 4.93 fold at 12 hpi, Table 3) as a pattern recognition receptor to sense viral dsRNA and triggers downstream antiviral reactions [65,66]. In addition, elevated expression of tripartite motif-containing protein 25 (TRIM25) at the very beginning of the infection event was believed to be related to the interaction with retinoic acid-inducible gene 1 (RIG-I) in detecting the viral RNA intermediates [67,68]. Apart from the regulation of RSAD2, the host cells membrane condition during caIBDV infection was probably influenced with cellular cholesterol metabolism. An over-expression of lipase A (LIPA) and cholesterol 25-hydroxylase (CH25H), together with the down-regulation of StAR-related lipid transfer protein 4 (STARD4), lanosterol synthase (LSS) and acetoacetyl-coA synthease (AACS) (Fig 6B), imply an evidence on the modulation of intracellular cholesterol contents in the infected cells. The synergetic effect of the regulation of these genes may alter the lipid raft arrangement and therefore inhibit membrane fusion event between the infected cells and virus, limit the clathrin-and cholesterol-dependent endocytosis hence prevent the propagation of the virus [69,70,71].

Potential strategies used by caIBDV to outwit host cell defenses
Viruses make use of various mechanisms to escape from the host defense actions. One of the viral strategies is to arrest the apoptosis initiation so as to facilitate viral replication process. Apoptosis is one of the host defense mechanism to minimize the spread of viruses. Our previous study demonstrated that apoptosis occurred at 48 hpi [21], while Jungmann et al. revealed that the IBDV-induced apoptosis was first observed at 12 hpi [72]. It has been well studied that NF-κB mediated apoptotic pathway is initiated upon viral infection. In uninfected cells, the dimer of NF-κB are sequestered in the cytoplasm by a family of κB inhibitors (IκB). Upon virus infection, these IκB proteins undergo signal-induced degradation by proteasome triggered by the activation of IκB kinase (IKK) and the NF-κB complex is released into the nucleus for the expression of specific genes leading into apoptosis and other antiviral functions. In this study, expression of both NF-κB1 and NF-κB2 were observed in mock-and caIBDV-infected DF-1 cells, but there was no significant changes between both groups at the time points tested (though the trends of both NF-κB proteins increased with time in caIBDV-infected group). Besides, the expression of two IKKs (IKK1 or CHUK and IKK2 or IKBKB) also showed no significant differences, implying that the degradation capacity of IκB between two groups were similar. The result, however, demonstrated that there was a distinct elevation of two IκB levels (NFKBIA or IκBα and NFKBIZ or IκBz) at 12 hpi in the caIBDV-infected group. The increase of these IκBs potentially hindered the function of IKKs and eventually reduced the free NF-κB amount, hence retarded the apoptotic process. Apart from the IκBs, the levels of tumor necrosis factor, alpha-induced protein 2 (TNFAI2) and 3 (TNFAIP3) were also provoked in the infected cells. It has been demonstrated that TNFAIP3 is a potent cellular inhibitor of NF-κB activation [73,74]. Furthermore, the caspases-mediated apoptosis was probably inhibited by the over-expressed baculoviral IAP repeat-containing protein 2 (BIRC2 or ITA) in the infected cells. ITA is a member of the inhibitor of apoptosis family that inhibit apoptosis by interfering with the activation of caspases [75,76]. Liu et al. suggested that the nonstructural protein (NS) of IBDV possessed anti-apoptotic function at the early stage of virus infection [77]. It is therefore reflecting that the NS of caIBDV regulates the expression of IκB, TNFAIPs and ITA in the infected cells which aims to delay the apoptosis in the first 12 hpi and reserve viable host cells for viral replication.
Switching of alternatively spliced EFR3A transcript isoforms RNA-Seq data reveals that two splice variants of EFR3 homolog A (S. cerevisiae) (EFR3A) transcript were identified in mock-and caIBDV-infected DF-1 cells respectively (Fig 4). It was shown that EFR3A-201 expressed in mock-infected cells, while the expression of EFR3A-202 isoform took over the previous one at 12 hpi in the caIBDV-infected cells. Amino acid alignment reveals that four residue differences located at the N-terminal immediately downstream to the first methionine residue between the two variants. EFR3A gene encodes a membrane protein [78], but the definite cellular function of EFR3A in chicken cells has not been reported yet. Odrowaz et al. reported EFR3A participated in the negative control of a ETS transcription factor ELK1 in a human epithelial cell MCF-10A [79]. Whereas it was shown to be involved to the epidermal growth factor receptor signaling pathway [80]. It was also suggested that the differential expression of EFR3A gene in auditory brainstem neurons of mice with hearing deficit [81]. The switch of splice variant expression may be involved in altering the lipid membrane condition, but the molecular mechanism of this deviations to host-virus interaction needs further characterization.

Conclusion
Interactions between IBDV and its host has been extensively studied. A number of researches revealed that infection of chicken cells, including fibroblast and bursal cells, with IBDV may lead into an elevation of cytokines and interferons at 1 dpi to 7 dpi periods in general. Virus-induced apoptosis via caspase-and NF-κB-mediated pathways were also demonstrated. This study, on the other hand, disclosed the early host-virus interactions. With the aid of RNA-Seq, a more comprehensive expression landscape was obtained. The result presents the events occurred before the elevation of downstream effectors. Apart from the regulators of cytokines and interferons, modulations were observed in the gene candidates involved in cell membrane fluidity. It is believed that the changes in membrane conditions contributes to the frontline host response against endocytosis of IBDV and hence prevent infections. On the other hand, the intensively expressed anti-apoptotic genes induced by caIBDV delay the programmed cell death and hence prolong the viral replication cycle in the host cells.
This study explore the initial host response in DF-1 cells upon caIBDV infection and the potential virus strategy to counteract with the host's action.