Analysis of the Early Immune Response to Infection by Infectious Bursal Disease Virus in Chickens Differing in Their Resistance to the Disease

ABSTRACT Chicken whole-genome gene expression arrays were used to analyze the host response to infection by infectious bursal disease virus (IBDV). Spleen and bursal tissue were examined from control and infected birds at 2, 3, and 4 days postinfection from two lines that differ in their resistance to IBDV infection. The host response was evaluated over this period, and differences between susceptible and resistant chicken lines were examined. Antiviral genes, including IFNA, IFNG, MX1, IFITM1, IFITM3, and IFITM5, were upregulated in response to infection. Evaluation of this gene expression data allowed us to predict several genes as candidates for involvement in resistance to IBDV. IMPORTANCE Infectious bursal disease (IBD) is of economic importance to the poultry industry and thus is also important for food security. Vaccines are available, but field strains of the virus are of increasing virulence. There is thus an urgent need to explore new control solutions, one of which would be to breed birds with greater resistance to IBD. This goal is perhaps uniquely achievable with poultry, of all farm animal species, since the genetics of 85% of the 60 billion chickens produced worldwide each year is under the control of essentially two breeding companies. In a comprehensive study, we attempt here to identify global transcriptomic differences in the target organ of the virus between chicken lines that differ in resistance and to predict candidate resistance genes.

I nfectious bursal disease virus (IBDV) is a highly contagious virus with a bisegmented double-stranded RNA (dsRNA) genome (belonging to the family Birnaviridae) which causes immunosuppression in chickens (1). Segment A contains two overlapping open reading frames, the larger of which encodes viral proteins VP2, VP3 (both structural capsid proteins), and VP4 (a viral protease). The smaller open reading frame encodes the nonstructural protein VP5. Segment B encodes VP1 which is a multifunctional polymerase. There are two known serotypes: serotype I viruses cause a range of disease severity in chickens and are further classified into classic, variant, and highly virulent strains, and serotype II viruses are nonpathogenic. Although largely controlled by vaccination, new virulent strains of the virus mean that infectious bursal disease (IBD; also known as "Gumboro" disease) still remain a threat to the poultry industry. The virus infects dividing IgM ϩ B lymphocytes and the main site of viral replication is the bursa of Fabricius, where B cells are produced (2,3). IBDV can also infect macrophages (3)(4)(5).
Infection is spread orally via contaminated feed and water (6). IBDV affects young birds, with the disease usually being diagnosed in 3-to 6-week-old birds. Younger birds do not show clinical signs but are immunosuppressed (7,8). Symptoms include anorexia, depression, diarrhea, ruffled feathers, immunosuppression, and bursal lesions. Death is often due to dehydration, which leads to kidney lesions (9). The disease peaks between 2 to 5 days postinfection (dpi) and is practically cleared by day 7 (10). During this acute phase, bursal follicles are depleted of B cells, and the bursa becomes atrophic. Abundant viral antigen can be detected in the bursal follicles and other lymphoid organs such as the cecal tonsils and spleen. CD4 ϩ and CD8 ϩ T cells accumulate around the site of virus replication (6). Mortality is variable and tends to affect layers more than broilers but can be up to 100% with very virulent strains of the disease. Even if birds survive, the resulting immunosuppression and effect on egg production in layer birds is significant (11).
Although disease is currently controlled by vaccination, alternative control measures are being explored, including breeding for enhanced genetic resistance. Involvement of the major histocompatibility complex (MHC) in resistance to IBDV has been a topic for debate, but it does appear to have a role (12). Two specific-pathogen-free lines of chickens held at The Pirbright Institute (previously known as the Institute for Animal Health), a Brown Leghorn line and the White Leghorn line 6 1 , were previously shown to be susceptible and resistant, respectively, to IBD (13). Differences in IBDV viral loads, measured by quantitative reverse transcription-PCR (qRT-PCR), in the bursae of infected birds of these two lines were evident from as early as 1 dpi (14). In the present study, we carry out a comprehensive gene expression study, comparing the response to infection in the two lines at 2, 3, and 4 dpi, with the hypothesis that genes underlying at least some of the resistance mechanisms will be involved at this early, innate stage of the immune response.

MATERIALS AND METHODS
Ethics statement. All animal work was conducted according to United Kingdom Home Office guidelines.
Experimental animals. Three-week-old birds from each line (6 1 and Brown Leghorn [Brl]) were separated into two experimental rooms, with ad libitum access to food and water. In one room, 54 birds were mock infected intranasally with 100 l of phosphate-buffered saline (PBS; 27 birds from each line). In the other room, 54 birds were infected intranasally (27 birds from each line) with 10 1.3 50% embryonic infectious doses (EID 50 ) of the classical virulent IBDV strain F52/70 in 100 l of PBS. The birds were monitored for clinical signs for up to 4 dpi and killed at 2, 3, or 4 dpi (nine individuals at each time point). Spleens and bursae were collected from all birds for qRT-PCR analysis for virus and host genes, for microarray analysis, and for bursal damage scores by immunohistochemistry. Blood samples were also taken from each bird for DNA isolation.
Bursal damage scoring. Segments of bursal tissue were fixed in formaldehyde-saline (pH 7.6) before routine histological processing and staining with hematoxylin and eosin. The degree of bursal damage was assessed using the histological scoring system (scored in the range 0 to 5, with 5 indicating the greatest level of damage) described by Muskett et al. (15), with each section being scored blindly by two separate individuals.
RNA preparation. Tissue samples (ϳ30 mg) were stabilized in RNAlater (Ambion/Life Technologies, Paisley, United Kingdom) and disrupted using a bead mill (Retsch MM 300; Retsch, Haan, Germany) at 20 Hz for 4 min. Total RNA was prepared using an RNeasy kit (Qiagen, Crawley, United Kingdom) extraction method according to the manufacturer's protocol. Samples were resuspended in a final volume of 50 l of RNase-free water. Concentrations of the samples were calculated by measuring the optical density at 260 nm (OD 260 ) and the OD 280 on a spectrophotometer (NanoDrop; Thermo Scientific, Paisley, United Kingdom). The quality of the RNA was checked on a bioanalyzer (Agilent Technologies, South Queensferry, United Kingdom). An RNA integrity number (RIN) of Ͼ8 proved the integrity of the RNA.
Microarray hybridization. Biotinylated fragmented cRNA was hybridized to the Affymetrix chicken genome array (Affymetrix, Santa Clara, CA). This array contains comprehensive coverage of 32,773 transcripts corresponding to Ͼ28,000 chicken genes. It also contains 689 probe sets for 684 transcripts from 17 avian viruses, including IBDV. For each experimental group (infected/control in two tissues at three time points in two lines), three biological replicates (three pools of RNA from three birds) were hybridized. Thus, 72 arrays were used in total. Hybridization was performed at 45°C for 16 h in a hybridization oven with constant rotation (60 rpm). The microarrays were then automatically washed and stained with streptavidin-phycoerythrin conjugate (SAPE; Invitrogen) in a Genechip fluidics station (Affymetrix). Fluorescence intensities were scanned with a GeneArray Scanner 3000 (Affymetrix). The scanned images were inspected and analyzed using established quality control measures. Array data have been submitted to Array Express (http://www .ebi.ac.uk/arrayexpress/) under the accession number E-TABM-1129.
Statistical analysis. Gene expression data generated from the GeneChip Operating Software (GCOS) were normalized using the PLIER (probe logarithmic intensity error) method (16) within the Affymetrix expression console software package. These normalized data were then analyzed using the limma and FARMS (17) packages within R in Biocon-ductor (18). Probes with a false discovery rate (FDR) of Ͻ0.05 and a fold change of Ն2 were deemed to be significant.
Analysis of differentially expressed genes. Gene ontogoly (GO) terms associated with genes differentially expressed (DE) during the host response were analyzed using EasyGO (http://bioinformatics.cau.edu.cn /easygo/). To determine which biological pathways are involved in the responses to viral infection, Pathway Express (http://vortex.cs.wayne.edu /projects.htm) was used. Genes differentially expressed during the host response (FDR Ͻ 0.05) were analyzed against a reference background consisting of all genes expressed in the experiment. Factors considered by Pathway Express include the magnitude of a gene's expression change and its position and interactions in any given pathway, thus including an "impact factor" when calculating statistically significant pathways. Anything with a P value of Ͻ0.25 is deemed significant when using this software.
Genes were clustered by similar expression pattern and analyzed for enriched GO terms and transcription factor binding sites (TFBS) using Expander (v5.2) (http://acgt.cs.tau.ac.il/expander/expander.html). Normalized expression data from control samples were compared to infected samples to examine the host response to IBDV infection. Enrichment analysis of particular GO terms or TFBS within clusters was done using the TANGO and PRIMA functions, respectively, within the Expander package. Use of the Ingenuity Pathway Analysis (IPA) program (Ingenuity Systems) revealed which canonical pathways and physiological functions were affected by IBDV infection in the host (Benjamini-Hochberg multiple testing correction; FDR Ͻ 0.05).
Viral quantification and specific gene expression analysis by quantitative real-time PCR. TaqMan real-time qRT-PCR was used to quantify viral RNA levels and for confirmation of the microarray results for the mRNA levels of selected genes. Primers (Sigma) and probe (PE Applied Biosystems, Warrington, United Kingdom) ( Table 1) were designed using Primer Express (PE Applied Biosystems). Briefly, the assays were performed using 2 l of total RNA and the TaqMan FAST Universal PCR master mix and one-step RT-PCR Mastermix reagents (PE Applied Biosystems) in a 10-l reaction. Amplification and detection of specific products were performed using the Applied Biosystems 7500 Fast real-time PCR system with the following cycle profile: one cycle at 48°C for 30 min and 95°C for 20 s, followed by 40 cycles at 95°C for 3 s and 60°C for 30 s. The data are expressed in terms of the cycle threshold (C T ) value, normalized for each sample using the C T value of 28S rRNA product for the same sample, as described previously (19)(20)(21). The final results are shown as 40-C T using the normalized value, or as fold change from uninfected controls.

Assessment of IBDV viral load in susceptible and resistant chicken lines.
After challenge of the two inbred lines, TaqMan analysis was used to measure viral load in bursal samples from control and infected birds from both lines. It became immediately obvious on analyzing the viral load data that the respective phenotype of the two lines (BrL susceptible and line 6 1 resistant) had reversed since the lines were last studied in previous decades (13,14). In repeated experiments, the line 6 1 birds were susceptible and the BrL birds were highly resistant, with only the odd BrL bird showing detectable viral load in the bursa postchallenge (Fig. 1A). In terms of bursal damage, all control birds had no signs of bursal damage, whereas all infected birds had some bursal damage ( Fig.  1B). At 3 and 4 dpi, high bursal viral loads in the line 6 1 birds correlated with high bursal damage scores, whereas bursal damage scores in infected BrL birds remained low (Fig. 2).
In order to confirm which biological processes were involved during infection of the birds and to determine whether there was any that had not been identified in previous, more-focused studies, we decided to analyze the gene ontology (GO) functional annotations of the genes being differentially expressed. The microarray data produced for the host response (infected versus control susceptible birds) in the spleen were analyzed with EasyGO (22), which examines DE genes for their association with particular GO terms compared to the microarray background as a whole. The response in the spleen was analyzed as opposed to that in the bursa since the large degree of cell damage occurring in the bursa would substantially mask the immune response. The genes differentially expressed during this experiment were seen to be involved in processes such as "immune response," "apoptosis," "response to stimulus," and "cytokine production," which is what would be expected in this type of infection.
The data were also analyzed using Pathway Express (23), which, based upon the KEGG pathways (24), pictorially illustrates the genes that are up-or downregulated in any given biological pathway. Figure 3 shows examples using the data from the spleen at 4 dpi. Genes involved in the extrinsic apoptosis pathway and the Toll-like receptor signaling pathway, which play integral roles during the innate immune response, are upregulated. In the bursa, genes involved in the B-cell receptor signaling and cell cycle pathways were dramatically downregulated (Fig. 4). Many of the biological responses seen in the bursa will be due to viral replication and cell damage and will not be antiviral responses per se. It must be borne in mind that these diagrams are based on the human pathways and so in some cases are not completely demonstrative of the chicken pathways, i.e., avian-specific genes are not represented. Other pathways seen to be significantly (FDR-corrected P value of Ͻ0.25) involved are shown in Table 3. In order to cluster the genes demonstrated to be involved in the host response to IBDV into groups with similar expression profiles, the CLICK algorithm within the Expander program (25) was used. The upregulated genes formed four distinct clusters, while the downregulated genes all cluster together in a fifth group.
Expander was then also used to look for enrichment of GO terms associated with genes, genomic locations and TFBS within the clustered DE genes (834 upregulated and 377 downregulated). The upregulated genes were over-represented in biological processes such as "immune response," "response to stimulus," "cytokine activity," and other functions associated with viral infection. Only "protein binding" was highlighted among the downregulated genes. Examination of the microarray data did not show any . IRF7 regulates transcription of type I IFN genes and IFN-stimulated genes by binding to ISRE motifs in their promoters. Genes containing a binding site for the Drosophila transcription factor Ovo (presumably the chicken orthologue thereof) was also over-represented. No TFBS were seen to be enriched within the downregulated array genes (Fig. 5).
Complementing the findings from Pathway Express, IPA was also used to identify the physiological functions and biological pathways most highly involved during the host response to IBDV infection. The most significant physiological functions are represented in Fig. 6A. It can be clearly seen that processes involving B and T cell development and differentiation constitute many of the pathways highlighted. Tissue development, necrosis, and mortality all seem to play an important role during IBDV infection. Pathways affected by IBDV infection are shown in Fig. 6B. Among the upregulated genes, those involved in signaling from several cytokine receptors (IFN, IL-6, IL-10, and GM-CSF/CSF-2), as well as in apoptosis, were differentially expressed. Genes involved in the activation of hepatic stellate cells were also differentially expressed. Interestingly, Ma et al. (26) recently described the effects of IBDV on Kupffer cells, macrophages found in the lining of the liver, suggesting that the liver is a site of IBDV infection and therefore possibly the innate immune response. Analysis of the downregulated genes highlights genes that are involved in endothelial cell development, proliferation, and migration.
Differences between susceptible and resistant lines. Resistance to IBD between the two lines could be due to a number of mechanisms. For example, BrL birds could simply express certain genes, such as those involved in key innate immune responses, at a constitutively higher level than line 6 1 birds, and thus mount a stronger innate response upon infection, limiting viral replication and thus disease. Alternatively, after infection, BrL birds could upregulate the expression of key immune function genes to a greater degree than line 6 1 birds, thus mounting a stronger induced immune response. Either or both mechanisms could contribute to IBD resistance.
There were genes that showed large expression differences between the two lines, even before infection, which in turn led to the differing responses seen upon infection. Examination of the control birds in the two lines showed considerable differences in expression levels of certain genes, including BLB1, SRFBP1, XRCC3, TNFRSF1B, GNG4, and IFITM3, all of which were more highly expressed in the resistant (BrL) than in the susceptible line (line 6 1 ). Some of these genes may therefore play an important role in disease resistance (382 DE probes in the spleen and 160 DE probes in the bursa: see Table S2 in the supplemental material).
Upon infection, differences in gene expression were also seen between the two lines. Genes more highly expressed in the resistant BrL line included BLB1, CARD9, BLVRA, AICDA, MYBL1, SFRP1, B6.3, VPREB3 and MMP13, whereas genes more highly expressed in the susceptible line included MMP9, IFNA, CCL5, MMP7, AVD, NOS2A, CXCLi2, IL-6, LYG2, GBP, CCL19, IF-ITM1, IFITM3, and IFITM5. Of course, "more highly expressed" could also mean "less downregulated." Again, it seems reasonable to assume that some of these genes may play an important role in resistance/susceptibility to IBDV. A full list of genes differentially expressed between the two lines in both the bursa (3,659 DE probes) and the spleen (1,570 DE probes) is given in Table S3 in the supplemental material.
Verification of the microarray results. Twenty-one genes were chosen for verification by qRT-PCR, including genes involved in the host response and genes differentially expressed between the susceptible and resistant lines (either inherently or during the course of infection) ( Table 4). Of the 20 genes tested, differential expression in the microarray was confirmed in the qRT-PCR analysis for 17. However, the microarray results for SRFBP1, TNFRSF1B, and CARD9 were not replicated in the qRT-PCR analysis (note that the qRT-PCR assay did not work for BLB1) (Fig. 7).
Potential candidate genes for IBDV resistance. In contrast to other avian pathogens, very little information is available regarding potential quantitative trait loci (QTL) regions for resistance to IBDV. The only published material is to be found in a Ph.D. thesis (14) that suggests that several genomic regions are potentially involved in resistance phenotypes, defined by death, bursal damage, IBDV genomic RNA levels, and IFN-␥ mRNA levels.
Due to the lack of robust accompanying genetic data, we had to identify potential candidate genes for resistance from our expression data alone. We found genes that are differentially expressed between susceptible and resistant lines either inherently or in response to IBDV infection. Examination of these genes, along with their known biology and how they might play a role in the pathways and biological networks identified in our analyses, allowed us to identify potential candidate genes for susceptibility/resistance to IBDV infection. Table 4 lists candidate genes (indicated in boldface) for resistance to IBDV, as determined by their differential expression between susceptible and resistant lines.

DISCUSSION
IBDV preferentially replicates in IgM ϩ B cells, which it enters after binding to cellular receptors, some of which have been identified, and induces host cell apoptosis, with both VP2 and VP5 playing a role in both binding and apoptosis. VP2 forms a subviral particle which binds to HSP90 (27). Other host molecules, including p53 binding protein (TP53BP1), stathmin (STMN1), and chondroitin sulfate, are also targets for VP2 (28). VP5 interacts with the volt- age-dependent anion channel 2 (VDAC2) protein (29) and with the p85␣ regulatory subunit of phosphatidylinositol 3-kinase (PI3K) (30). Integrin ␣4␤1 is also a specific host binding receptor for IBDV (31). Comparing gene expression levels between control birds of each line, we found higher mRNA expression levels of some of these receptor genes (HSPCB [HSP90B], TP53BP1, STMN1, and ITGB1) in the resistant line (BrL) than in the susceptible line (6 1 ). It seems counterintuitive that the resistant line should have more receptors for the virus than the susceptible line. This presumably reflects some mechanism whereby viral entry into the cell is limited. The Jun NH 2 -terminal kinase and p38 mitogen-activated kinase pathways (30,32) have previously been suggested to be involved in the pathogenesis of IBD. We confirm here the involvement of these pathways during IBDV infection and highlight the whole spectrum of genes in these pathways whose expression is altered. The importance of apoptosis in IBD pathogenesis is clear from the EasyGO analysis and Pathway Express highlights the affected molecules, also showing upregulation of PI3K. IBDV activates PI3K/Akt signaling through binding of the nonstructural VP5 protein to the p85␣ regulatory subunit of PI3K, resulting in the suppression of premature apoptosis and improved virus growth after infection (30). In the bursa, there is also considerable downregulation of genes involved in the B cell receptor signaling pathway and in the cell cycle. The former probably reflects the huge reduction in B cell numbers as they are destroyed by the infection. It remains to be seen whether VP5, or indeed any of the other IBDV proteins, also targets the cell cycle. Generally, T cells are refractory to IBDV infection with IBDV but promote virus clearance, by mechanisms that are not well understood. Rauf et al. (33) showed that CD4 ϩ and CD8 ϩ T cells enter the infected bursa and that cytotoxic T cells play a role in clearing infected cells through the perforin-granzyme pathway, as identified by the upregulation of mRNA expression levels for perforin and granzyme, DNA repair and apoptotic proteins, high-mobility protein group and poly(ADP-ribose) polymerase (PARP) in the bursa, and the presence of perforin-and granzyme-expressing CD8 ϩ T cells. Similarly in the present study, we observed upregulation of granzyme A, various PARP genes, and also IFN-␥. Rauf et al. (33) also noted a decrease in NK-lysin mRNA expression levels, suggesting a reduced role for NK cells. In contrast, the present study demonstrated a 5-fold upregulation of NK lysin in the bursa at 3 dpi, rising to 15-fold at 4 dpi, suggesting that in these birds both cytotoxic T cells and NK cells are involved in the response to the virus. The chicken has a somewhat different repertoire of pattern recognition receptors (PRR), including Toll-like receptors (TLR), to that of mammals (34), but the repertoire of potential pathogenassociated molecular patterns (PAMP) recognized by these PRR is thought to be similar. Little is understood of the effects infection with different classes of pathogens have on the expression repertoire of the different TLR. In the present study, during the early stages of IBDV infection in the bursa, mRNA expression levels of TLR1LB, TLR2A, TLR3, TLR4, and TLR15 (and TLR21 in the spleen) and the cytosolic PRR MDA5 were upregulated, whereas those of TLR1LA, TLR2B, and TLR7 were downregulated. Other studies have also reported this differential regulation of TLR3 and TLR7 following infection with a classical strain of IBDV (33,35). TLR3 recognizes viral dsRNA, whereas TLR7 recognizes viral single-stranded RNA; the genome of IBDV is dsRNA, and upregulation of TLR3 might represent a positive feedback loop to drive the innate inflammatory antiviral response. MDA5 also recognizes viral dsRNA, and its mRNA expression levels were also upregulated (70-fold at 3 dpi compared to 9-fold for TLR3) in response to IBDV infection.
It is more difficult to explain the differential mRNA expression patterns of some of the TLRs that are thought to generally recognize components of the surface of pathogens, particularly bacteria. The precise ligand of the chicken-unique TLR15 remains to be determined. It has been described to recognize a variety of PAMP, including diacylated lipopeptides (36) and a yeast-derived agonist (37). We previously reported increased expression of TLR15 after infection with Marek's disease virus (38). Together with the data presented in the present study, this implies a role for TLR15 in antiviral responses and that it perhaps recognizes a surface component of viruses. It is interesting to note the differential response of the TLR1 and TLR2 paralogues, i.e., TLR1LB and TLR2A mRNA expression levels are upregulated, and those of TLR2B and TLR1LA are downregulated. Moreover, TLR2-1 and TLR1-2 heterodimers cooperatively signal the presence of PGN, diacylated lipopeptides and MALP-2, whereas TLR2-1 and TLR1-1 heterodimers did not recognize MALP-2. Both combinations, however, recognized the triacylated lipopeptide, Pam3 (39). It is therefore clear that the two different heterodimer combinations discriminate between different ligands and that they might differentially recognize surface components of viruses. Analysis of the samples from IBD-susceptible and -resistant birds, with or without infection with IBDV, using whole-genome microarrays identified genes that either show different inherent levels of gene expression (without infection) between the lines or that are transcribed differently after infection, thus potentially eliciting differing host responses. They are thus prime candidates for future testing in genetic mapping studies, either as targets for knockout experiments or as targets for direct interaction with IBDV proteins.