Transcription profiles of the responses of chicken bursae of Fabricius to IBDV in different timing phases

Infectious bursal disease virus (IBDV) infection causes immunosuppression in chickens and increases their susceptibility to secondary infections. To explore the interaction between host and IBDV, RNA-Seq was applied to analyse the transcriptional profiles of the responses of chickens’ bursas of Fabricius in the early stage of IBDV infection. The results displayed that a total of 15546 genes were identified in the chicken bursa libraries. Among the annotated genes, there were 2006 and 4668 differentially expressed genes in the infection group compared with the mock group on day 1 and day 3 post inoculation (1 and 3 dpi), respectively. Moreover, there were 676 common up-regulated and 83 common down-regulated genes in the bursae taken from the chickens infected with IBDV on both 1 and 3 dpi. Meanwhile, there were also some characteristic differentially expressed genes on 1 and 3 dpi. On day 1 after inoculation with IBDV, host responses mainly displayed immune response processes, while metabolic pathways played an important role on day three post infection. Six genes were confirmed by quantitative reverse transcription-PCR. In conclusion, the differential gene expression profile demonstrated with RNA-Seq might offer a better understanding of the molecular interactions between host and IBDV during the early stage of infection.


Background
Infectious bursal disease virus (IBDV) is an important small RNA virus in the family of Birnaviridae, which may lead to severe immunosuppressive effects and pathological damage in young chickens. It mainly targets early B cells, especially those in the gut-associated lymphoid organ, the bursa of Fabricius [1]. After infection with IBDV, the bursa collected from chickens displays oedema and haemorrhage, or even necrosis, accompanied by a large volume of T cells infiltration and a cytokine storm [2]. Especially when chickens mix infection with IBDV and other pathogens, such as E. coil and Newcastle disease virus, farms can suffer large economic losses. Therefore, it is urgent to learn the pathways of IBDV invasion and replication in host cells and the interactions between IBDV and host, in order to determine effective strategies for the prevention and control of IBDV infection.
Transcriptional profiling analysis has displayed a great deal of information about host-pathogen interactions, including many important biological processes, cellular components and molecular functions. Thomas Ruby et al. have analysed gene expression levels in the bursae of young chickens from the resistant and susceptible lines by using a novel cDNA microarray and shown that the changes in gene expression in the target tissue during the early stages of infection of young chickens with IBDV [1]. Similarly, Wong RT et al. have screened differentially expressed genes in IBDV-induced apoptotic chicken embryonic fibroblasts by using cDNA microarrays and unravelled the candidate physiological pathways involved in host-virus interactions on a molecular level in vitro [3]. However, cDNA microarray is a closed system which only detects those specific genes. It can not offer comprehensive transcriptional landscape described in the cells and bursae post IBDV infection. We still do not entirely understand how the host defends against IBDV and how the infection influences the biological metabolism of the host at different stages of infection. With the development of high throughput sequencing, it is feasible to obtain transcript splice-variants, isoformas and new genes and compare the biological metabolism and immune responses of hosts at different infection points [4], which will provide new insights for understanding the pathogenesis of IBDV and the antiviral immunity of the host.

Birds and viruses
A classical IBDV strain CJ801 was kindly provided by Prof. Jue Liu of the Beijing Academy of Agriculture and Forestry, Beijing, China. All the studies were performed on 3-week-old specific-pathogen-free (SPF) chickens (Vital Bridge Co. Ltd., Beijing, China). The experimental procedures and animal management were approved by the Institutional Animal Care and Use Committee at Henan Institute of Science and Technology.

Experiment design
Eighteen SPF white leghorn chickens were randomly divided into two groups with 9 chickens for each group: the mock group (the healthy group) and the IBDVinoculated group (the infection group). Chickens from the infection group were inoculated with IBDV CJ801 stock at a dosage of 10 3 EID 50 /0.1 mL through eye-nose drops [5,6]. The chickens of the mock group were mock challenged with PBS. All chickens were kept in two separate isolator and they were free to drink and eat with 12 h light and 12 h dark. On days 1, 3 and 7 after infection, 3 chickens of each group were randomly selected and killed under anesthesia for bursa collection [7]. Each bursa was immediately put into liquid nitrogen and then stored in 80°C refrigerator. RNA sequencing was performed with total RNA from bursae from each group at the first two time points and completed by a commercial company (Shanghai Biotechnology Corporation, Shanghai, China).

RNA extraction and purification
Total RNA was separately extracted from each bursa using RNAiso Plus Total RNA extraction reagent (TAKARA, China) following the manufacturer's instructions and checked for RNA integrity number (RIN) to inspect RNA integrity by an Agilent 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, US). Qualified total RNA was further purified by RNeasy Micro Kit (QIAGEN, GmBH, Germany) and RNase-Free DNase Set (QIAGEN, GmBH, Germany). The qualified RNA meets the following requirements: the RIN value was not less than 7.0, while the ratio of 28 s/18 s should be between 1.8 and 2.2.

cDNA library construction and sequencing
We used at least 3 μg of high-quality mRNA for each library construction. The mRNA was fragmented by incubating at 94°C in fragmentation buffer to yield a size range of 400-500 bp, verified by an Agilent 2100 Bioanalyzer (Agilent). We used SuperScript II Reverse Transcriptase (Invitrogen) for the library construction. According to the corresponding process shown in the cBot User Guide, we completed the Cluster generation and the first sequencing primer hybridisation on cBot in Illumina HiSeq 2500 sequencer. Then, sequencing was performed using the Illumina HiSeq 2500 according to its standard protocols. The data for each sample was not less than 4.0 G while the percentage of Q20 (bases of Q > =20/all bases of sequencing) was not less than 90%.

Data analysis
After the sequencing was completed, the obtained image data was transformed into raw reads and stored in a FASTQ format. The raw reads were cleaned by a shortreads pre-processing tool (FASTX-Toolkit, version 0.0.13). The low-quality reads, including adapter, ribosome RNA and sequences shorter than 25 nt with Q < 20 at 3′ end, were removed. The resultant clean reads from each sample library were used for further analyses. The clean reads were mapped to the chicken genome (galGal4, Ensembl release 85) by a mapping tool TopHat2 (version: 2.0.9). The differentially expressed genes with fold changes ≥ 2 or ≤ 0.5 and False Discovery Rate (FDR) < 0.05 were analyzed using the web-based tools in DAVID to identify enriched gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, group functionally related genes and cluster the annotation terms. The raw sequencing data have been deposited in the Gene Expression Omnibus (GEO) at the NCBI under accession number GSE94500.

Confirmation of RNA-Seq data by quantitative real-time PCR
To verify the correctness of the RNA-Seq results, some genes were randomly selected for quantitative real-time PCR (qRT-PCR) experiments. CCR2 and CCR4 were used for their chemotaxis effects on monocyte, macrophage and dendritic cells. IL-6 is an important inflammatory factor. IFN-γ is over-expressed during IBDV infection and plays an important role in host antiviral responses. MHC-I molecule is to display intracellular protein to cytotoxic T cells. Like IL-2, IL-15 exerts a huge role in the process of virus infection. So, CCR2, CCR4, IL-6, IFN-γ, MHC-I, IL-15 were chosen in the study while β-actin was used as housekeeping gene [8].
The primers used in this study are listed in Table 1. Total RNA was extracted from chicken bursae by total RNA Kit (Omega Bio-Tek, Guangzhou, China) and then reverse transcribed into cDNA using M-MLV Reverse Transcriptase (Promega, CA, USA) following the manufacturer's instructions. The qRT-PCR was performed using a 7300 Real-Time PCR System (ABI, Britain) with SYBR® Select Master Mix (ThermoFisher Scientific, China). The reaction conditions were 30 s at 95°C, 40 cycles of 5 s at 95°C, 25 s at 58°C and 30 s at 72°C, followed by a final step of melting curve analysis. After amplification, the relative fold change of the differentially expressed genes was calculated through the 2 −ΔΔCT method [9]. All samples were performed in 3 replicates as technical repeats to ensure the reproducibility of the amplification.

VP2 gene expression during IBDV infection in chickens
After being challenged with IBDV, chickens from the infected group displayed some clinical symptoms on 1 dpi, such as decreased feed intake, fluffed feathers and diarrhoea. Autopsy results displayed haemorrhaging and necrosis on the bursae of Fabricius. However, chickens from the mock group displayed very healthy. To further confirm the success of viral infection, IBDV VP2 genes in the bursa from the infected group were quantified by real-time PCR on 1 and 3 dpi. The results showed that VP2 gene expression significantly increased in the infected group and was about 21 and 79 times higher than those of the mock group on 1 and 3 dpi, respectively ( Fig. 1).
Identification of differentially expressed genes between the mock group and the infection group on day 1 and day 3 To analyse the gene expression of IBDV-infected chickens, the total RNA was prepared from the bursae of Fabricius of IBDV or mock-infected chickens on days 1 or 3 after infection. In this study, we identified a total of 15546 genes in the chicken bursa libraries. As shown in Fig. 2a, among the annotated genes, we identified 2006 and 4668 differentially expressed genes in the infection group compared with the mock group on 1 and 3 dpi, respectively (FDR ≤ 0.05). Out of these, 1618 genes were the common differentially expressed genes with 388 and 3050 specific differentially expressed genes, respectively, on 1 and 3 dpi. Among the differentially expressed genes, there were 946 up-regulated and 179 down-regulated significantly differentially expressed (SDE) genes on 1 dpi while 2142 up-regulated and 1597 down-regulated SDE genes were identified on 3 dpi in the infection group with a foldchange ≥ 2 or ≤ 0.5 ( Fig. 2b-c). Among these SDE genes, there were 676 common up-regulated and 83 common down-regulated genes in the bursae taken from the chickens infected with IBDV on both 1 and 3 dpi (Fig. 2b-c).

The host mainly displayed strong immune responses in the early stage of virus infection
The functions and pathways of the 799 differentially expressed transcripts were analysed based on the Gene Ontology (GO) project and KEGG using the SBC Analysis system. Top gene ontology clusters of the common SDE genes between the infection group and the mock group were mainly grouped into molecular function and biological processes, such as binding, response to stress, external stimulus, biotic stimulus, external biotic stimulus, virus, lipopolysaccharides and molecules of bacterial origin, immune responses and cytokine activity ( Table 2).
KEGG pathway enrichment analysis showed that the differentially expressed genes were mainly clustered into the following pathways: cytokine-cytokine receptor interaction, influenza A, Herpes simplex infection, cell adhesion molecules and some signaling pathways (such as JAK-STAT, Toll-like receptor and RIG-I-like receptor) ( Table 3).

Validation of the RNA-Seq data by qRT-PCR
To validate the results of the RNA sequencing, qRT-PCR was carried out to determine the differentially expressed genes. Six (CCR2, CCR4, IL-6, IFN-γ, MHC-I, IL-15) genes were chosen for the qRT-PCR tests. The qRT-PCR results showed a similar expression pattern to that observed in the RNA sequencing analysis, even though the fold changes measured by the two methods were not entirely the same (data not shown, see Additional file 1: Table S2). The results revealed that the RNA sequencing results could represent all the gene expression variations.

Discussions
IBDV can cause immune suppression of poultry and bring great economic losses for the poultry industry. The host generally displays a series of antiviral responses to IBDV infection, including natural and acquired immunity. Moreover, the host's responses have different characteristics in different stages of infection. Therefore, studying the hosts' immune responses in different stages of infection is very important for better understanding of the pathogenesis of pathogenic microbes. It is both reasonable and feasible to obtain these differentially expressed genes by the RNA sequencing along with the progress of transcriptome sequencing technology. In this manuscript, we report for the first time the transcriptome changes of chicken bursae during the early stage of IBDV infection by RNA-seq. The results displayed that, even if all samples were taken from the infection group, gene changes were not exactly the same for the days 1 and 3 post inoculation. Meanwhile, there were differentially expressed genes between the infection group and the mock group post-infection.
Following infection and replication of IBDV, lots of T cells infiltrate the bursae of infected chickens [10], while IBDV causes injury to B cells and macrophages [11], thus inducing a so-called "cytokine storm", such as proinflammatory cytokines, antiinflammatory cytokines, chemokines, interleukins, nitric oxide (NO) and so on. Obviously, it was no doubt that the primarily observing immune responses by host were displayed by the primary or secondary infected cells on 1 dpi. On 3 dpi, these immune signals might still be present, but they were being diluted and confounded by the responses of different (myeloid & lymphoid) cell types, some of which will be displaying their own activated phenotype. Therefore, it was difficult to work out the signals coming from the resident stromal or immigrating immune cells, especially with regard to the activation of those cells. In this study, though similar results were shown in IBDVinfected DF1 cells at the early stage of infection [12], significant elevations of some important genes expression levels were displayed in the infection group and inflammatory response genes, antiviral related genes and these characteristic differentially expressed genes on 1 and 3 dpi were mainly observed.

Inflammatory cytokine storm during the stage of IBDV infection
Inflammatory response genes, such as NOS2, IL6 and TNFRSF1B, were up-regulated while some positive regulation of inflammatory response genes, such as TLR3, STAT5 and LRRC32, were also elevated post-infection with IBDV. The chemokines (CCL4, CCL19, CXCL12) and their receptors (CCR5, CCR7 and CXCR4) were involved in inflammatory responses. Among these genes, LRRC32 (leucine rich repeat containing 32, also known as GARP) is critical for tethering TGF-β to the cell surface, and rearrangements of LRRC32 or a neighbouring gene may be important for the pathogenesis of hibernomas [13,14]. The main sources of these chemokines are dendritic cells, T cells or white cells and their effector cells are T cells or other immune cells [15,16]. Thus, these genes might be involved in virus antigen presentation and virus killing during IBDV infection.
Tumour necrosis factor α (TNF-α) is a pleiotropic cytokine and could respond to a wide range of stimuli [17]. Along the TNF-mediated signalling pathway, several gene expressions were largely influenced during IBDV infection, such as TNFSF10 and TNFAIP2. In the present study, TNF receptor-associated factor 2 (TRAF2) was differentially up-regulated by two times compared to the mock group and has been shown to be essential for TNF-α-mediated activation of JNK. This could contribute to the activation of NF-kappaB and antiapoptotic signals [18]. Moreover, tyrosine phosphorylation of signal transducer and activator of transcription 2 (STAT2) have been reported to be induced by activation of JAK kinases because of type I IFN binding to cell surface receptors, leading to activating the expression of interferon-stimulated genes and driving the cell into an antiviral state [19]; it was also increased during IBDV infection. Tumour necrosis factor receptor superfamily member 4 (TNFRSF4), a receptor for TNFSF4/OX40L/ GP34, is a costimulatory molecule implicated in longterm T cell immunity and could act as a receptor for human herpesvirus 6B [20,21]. Tumor necrosis factor alpha-induced protein 2 (TNFAIP2) functions as an important pro-inflammatory gene and plays vital role in the process of inflammatory response [22].

Strong host antiviral immune responses were activated by IBDV infection
After IBDV infection, some anti-infection genes were significantly changed without regard to the first day or third day post inoculation, including defense response to virus, positive regulation of T cell-mediated cytotoxicity   Table 3 Top KEGG pathways associated with SDE genes between the mock group and the infection group  promoter regulator of the IFN-α/β genes and MHC-I antigens in chicken fibroblast cell line C32 [24,25], IRF-1 was up-regulated to 12-and 7-fold changes in the infection group on days 1 and 3, respectively, while the fold change of the IRF-1 gene expression in CEFs was 23-fold in a previous study on day 3 following infection with IBDV [26]. GATA binding protein 3 (GATA3) belongs to the GATA family of transcription factors and plays an important role in T cell development. It has been displayed to improve the secretion of IL-4, IL-5 and IL-13 from Th2 cells and induce Th0 cells differentiation into a T cell subtype [27]. Therefore, overexpression of GATA3 might be involved with virus inhibition of infiltrated T cells in this study. SAM domain and HD domain-containing protein 1 (SAMHD1), is a cellular enzyme in dendritic cells, macrophages and monocytes, responsible for blocking the replication of HIV by depleting the intracellular pool of deoxynucleoside triphosphates [28][29][30][31]. Radical S-adenosyl methionine domain-containing protein 2 (RSAD2, also known as CIG5, viperin or CIG33) with 361 amino acids is usually upregulated and displays antiviral defence responses against some pathogens, such as bovine respiratory syncytial virus (BRSV) and hepatitis C virus [32,33].
Host responses displayed significant differences on days 1 and 3 post inoculation In our study, on day 1 post-infection, the cytokinecytokine receptor interaction pathway, JAK-STAT signaling pathway, AGE-RAGE signaling pathway in diabetic complications, Toll-like receptor signaling pathway and intestinal immune network for IgA production were activated and related genes were largely up-regulated in the infection group (Fig. 3). However, on day 3 postinfection, many metabolism pathways were activated and participated in the process of viral infection. Wong has described that transcription profiles of IBDV-infected cells were also different at 24, 48, 72 and 96 h postinfection, which was in accordance with our results [3]. Hui et al. also reported that different isoforms were involved in IBDV infection at different time-points postinfection, such as the IFIT5-IRF1/3-RSAD5 pathway in the DF1 cells that was triggered in the early infection stage [12]. Early host responses were vital for prevention of viral infection and replication. Our study showed that CSF3 (Myelomonocytic growth factor) increased in the bursae of IBDV-infected chickens at 1 day post infection, which could induce leukotrienes and exert anti-HIV-1 effects [34]. FOS (FBJ murine osteosarcoma viral oncogene homolog), as a cellular immediate-early gene was markedly induced before expressions of the Epstein-Barr Virus (EBV) transactivator genes were activated [35], which was similar in the process of IBDV infection. Likewise, host immune-related genes, such as VCAM1, NFKB1 and TLR2/3, were also elevated at the early stage in the infection group, which were consistent with many previous reports [36][37][38].
In the current study, several metabolic pathwaysassociated genes were involved on day 3 post-infection with IBDV, such as MYL9, GSTA and GMPPA, and these corresponding genes changed significantly compared with those on day 1. Among them, many metabolic pathways played vital roles in tissue injury and cytoskeleton repair, which might meet higher metabolic demands for virus replication and bursal injury repair during the middle stage of IBDV infection. It has been suggested that glycosaminoglycan could participate in a variety of biological processes, including cell-matrix interactions and activation of chemokines, enzymes and growth factors, and be used as a candidate to enhance chronic wound repair [39,40]. Sun et al. have displayed that an array of factors involved in glycerophospholipid metabolism, arachidonic acid metabolism and tryptophan metabolism pathways were activated during human hepatitis C virus infection and arachidonic acid metabolism was also an up-regulated marker of neuro-inflammation in HIV-1 transgenic rats [41,42]. Moreover, glycosphingolipid biosynthesis played an important role in nephritis development and CD1d-mediated lipid antigen presentation [43,44]. Glutamine is used in the tricarboxylic acid cycle and human cytomegalovirus infection could activate the mechanisms that switch the anaplerotic substrate from glucose to glutamine to meet the biosynthetic and energetic needs of the viral infection [45].

Conclusions
Host (including fibroblast, DF1 and bursal cells) immune responses against IBDV have been widely reported, and many genes involved in viral infection have been studied. However, there were few reports on bursae-virus interactions in the early stage of IBDV infection by RNA-Seq. The results presented in this study have shown that lots of T cell infiltrated into bursae post IBDV infection accompanied with host immune responses and inflammatory responses. Overall, the host mainly displayed antiviral responses in the early stage of viral infection while many metabolic pathways were involved in the middle stage of viral infection. Therefore, finding an effective method of promoting IBDV-infected cells apoptosis and improving antiviral responses might be the best way to delay virus replication in the early stage. Moreover, inhibition of metabolic pathways may be an effective method of inhibiting the spread of the virus and the formation of a large amount of connective tissue in the middle or latent stages of viral infection.