Bursal transcriptome profiling of different inbred chicken lines reveals key differentially expressed genes at 3 days post-infection with very virulent infectious bursal disease virus

Infectious bursal disease is a highly contagious disease in the poultry industry and causes immunosuppression in chickens. Genome-wide regulations of immune response genes of inbred chickens with different genetic backgrounds, following very virulent infectious bursal disease virus (vvIBDV) infection are poorly characterized. Therefore, this study aims to analyse the bursal tissue transcriptome of six inbred chicken lines 6, 7, 15, N, O and P following infection with vvIBDV strain UK661 using strand-specific next-generation sequencing, by highlighting important genes and pathways involved in the infected chicken during peak infection at 3 days post-infection. All infected chickens succumbed to the infection without major variations among the different lines. However, based on the viral loads and bursal lesion scoring, lines P and 6 can be considered as the most susceptible lines, while lines 15 and N were regarded as the least affected lines. Transcriptome profiling of the bursa identified 4588 genes to be differentially expressed, with 2985 upregulated and 1642 downregulated genes, in which these genes were commonly or uniquely detected in all or several infected lines. Genes that were upregulated are primarily pro-inflammatory cytokines, chemokines and IFN-related. Various genes that are associated with B-cell functions and genes related to apoptosis were downregulated, together with the genes involved in p53 signalling. In conclusion, bursal transcriptome profiles of different inbred lines showed differential expressions of pro-inflammatory cytokines and chemokines, Th1 cytokines, JAK-STAT signalling genes, MAPK signalling genes, and their related pathways following vvIBDV


INTRODUCTION
Infectious bursal disease virus (IBDV) is a non-enveloped, bi-segmented dsRNA virus, a causative agent of infectious bursal disease [IBD (Gumboro disease)], an economically important immunosuppressive disease affecting young chickens between 3-6 weeks old [1].The virus infects and destroys primarily B-cells found in the bursa of Fabricius (BF) and macrophages [2].A recent study also showed that the virus infects and destroys dendritic cells [3].IBDV can be classified as vaccine, classical virulent and very virulent (vv) strains based on the virus pathogenicity [4].Although the first outbreak of vvIBDV was reported in Europe in the early 1990s [5], vvIBDV has been widespread and detected in the majority of poultry farms in Asia, Africa and South America [6].Clinical signs induced by the vv strain were similar to those induced by the classical virulent strain.However, the vvIBDV induced IBD with a more Downloaded from www.microbiologyresearch.orgby IP: 54.70.40.11On: Sun, 30 Dec 2018 08:09:22 pronounced bursal lesion and higher mortality rate compared to other IBDV strains [7].In addition, vaccination against the vv strain encounters many challenges as vvIBDV can break protective antibody titres, and therefore requires more efficient vaccination approaches [8].
Chickens through the development of diverse inbred lines have been used to investigate the genetic response mechanisms to different viral diseases such as Marek's disease (MD), avian leukosis and IBD [9][10][11][12].Specific-pathogenfree (SPF) layer-type chickens have been shown to be more susceptible to both classical and vv strains of IBDV compared to broiler-type chickens based on the clinical disease, mortality, bursal lesion score and antigen load [7].In addition, different inbred chicken lines exhibit different susceptibility levels based on the mortality rate, with lines 15 and 7 considered as highly resistant (0 % mortality), lines N and 6 as less resistant (8 % mortality) and line BrL as highly susceptible (80 % mortality) following infection with vvIBDV strain CS89 [12].Microarrays, in particular, have been widely used in characterizing the immunogenomic response against avian viral diseases, namely MD, IBD, avian influenza and Newcastle disease [13].Microarray analysis of susceptible line BrL and resistant line 6 following infection with F52/70 IBDV strain demonstrated early changes in the gene expression of the bursa, which may influence the mechanisms that regulate the differential genetic resistance [14].However, a subsequent study using a similar virus strain reported a contradicting result with line BrL being resistant and line 6 being susceptible based on the viral load and bursal damage, with the resistant line having lower viral load and lower bursal damage compared to the susceptible line [11].The discrepancy between these two studies is not clear but is probably due to the timing differences, whereby the first study was carried out a decade ago.Nevertheless, both studies emphasized the importance of an early innate immune response for the resistance mechanism [11,14].However, the transcriptome analysis of bursal tissue during peak IBDV infection is lacking.A previous study has shown that Rhode Island Red (RIR) chickens are susceptible to vvIBDV infection [15]; however, no study has been carried out on different inbred chicken lines.Only recently, next-generation sequencing (NGS)-based transcriptome profiles of SPF chickens following infection with classical IBDV were characterized [16].
In this study, genome-wide host responses following vvIBDV infection in bursal tissues of inbred chicken lines, namely lines 15, 6, 7, N, O, and P, were characterized based on the differential expression profiling using the strandspecific RNA-sequencing (RNA-seq) approach and available chicken genome sequences.The aim of this study is to characterize the bursal transcriptome of different inbred chicken lines following vvIBDV infection at 3 days postinfection (p.i.), in the effort to identify genes and their interactions in modulating the severity of the vvIBDV infection.This study provides new insight into the host-pathogen interaction in different inbred lines, highlighting important genes and pathways, and identifies potential biological modulators to be utilized as therapeutic agents for future vaccine development strategy against IBDV infection.

IBDV infection
The birds were monitored and evaluated for clinical disease development based on the clinical scoring system in Table S1 (available in the online version of this article).Throughout the course of the experiment, none of the control birds showed clinical signs of the infection.On the other hand, the clinical signs started to appear at 2 days p.i. and worsen at 3 days p.i. in the infected groups.The birds require euthanization when the total clinical score reached (+5) for 2 consecutive days.By 3 days p.i., all infected birds were euthanized since they showed clinical scores of (+5) for 2 consecutive days.The severity of IBDV infection across six different inbred lines was further evaluated based on the histological damage of the BF.The bursa of all infected birds showed severe bursal damage ranging between 3.5-4.4,with line 6 having the highest average bursal score and line O the lowest average bursal score (Table 1).There was no evidence of bursal lesion in the mock-infected bursa.

RNA-seq for cDNA and External RNA Controls Consortium (ERCC) libraries
The cDNA libraries generated by Illumina HiSeq 2000 successfully produced more than 52 million reads for each of the chicken lines with a library size between 150-200 nucleotides.Approximately, 50 million clean reads were generated with 96 % of them showing read lengths of 70-79 nucleotides.The ERCC spike-in RNA transcripts quantified were plotted as a dose-response curve and R 2 was determined from known ERCC transcript number in relation to the read density output (fragments per kilobase of transcript per million mapped reads, FPKM).These results indicate linear regression with R 2 obtained across all chicken lines is approximately 0.93 on log 2 -transformed plot, showing a strong correlation between the transcriptome read counts and RNA inputs.In addition, the R 2 between the replicates of each line are very similar with maximum differences of only 0.02 (Table S2).The detection limits for ERCC dose response is 0.000143 attomole µg À1 of the total RNA depicting the sensitivity of the NGS platform, and as recommended by the manufacturer, 2-5 % of the reads generated were successfully mapped to the ERCC reference sequences.
ERCC fold-change response assessment results reported the fold-change distribution of six inbred lines with R 2 between the range 0.75-0.85.The sequencing background noise detected from the dose-response curve was between the range 1.87E-5-8.13E-5FPKM, indicating that the platform used in this study has low background noise (Table S2).Generally, more than 14 000 transcripts were expressed in all lines.Details on the transcriptome profiling and the read mapping statistics are presented in Table S3.

Viral load quantification and read mapping to the IBDV genome
Viral load quantification was determined in each sample using real-time quantitative PCR, and the results showed that lines 6, 7, O and P have slightly higher viral load of 1.3fold compared to lines N and 15 (P<0.05,Fig. 1).Sequence reads that were mapped neither to chicken reference genome nor ERCC spike-in transcripts were retrieved and mapped to the IBDV genome (NC_004178 and NC_004179) [17].The mapping results showed that the mock-infected samples did not contain sequences that correspond to the IBDV genome, indicating they were not exposed to IBDV, whereas an average of 20 000-66 000 reads (0.6-2.3 % of reads used for mapping) from infected samples were successfully mapped to the IBDV genome (Table S4).However, the percentage of reads mapped to the viral genome did not correlate with the viral load quantified, suggesting that the RNA-seq data could not be used to infer viral load in an infected host.

Identification of differentially expressed (DE) genes
The correlation analyses of the sample replicates were performed for each chicken line to verify its reproducibility prior to differential expression analysis.The lowest correlation detected was in line O with R 2 =0.83, while the highest correlation was detected in line P with R 2 =0.95 based on log 10 expressed transcripts (Fig. 2).Gene expression of annotated chicken genes were measured from a total of 17 936 reference genes in the NCBI database, and ~80 % of the reference genes were successfully detected with a cutoff value of FPKM !1E-5.The distribution of the expressed genes is depicted as a volcano plot (Fig. S1).
A false discovery rate was used to control multiple analysis and also to determine the threshold of the Q-value, an analogue of the P-value [18].The DE genes were defined as genes that were measured to be at log 2 fold-change <À2 or log 2 fold-change >2 between the infected and mock-infected samples with Q-value <0.05.In addition, genes with expression detected only in one condition; either mock-infected or  infected samples with Q-value <0.05 were considered as uniquely expressed.Under this rule, 4588 genes were successfully identified to be DE genes (Table S5).The summary of DE genes in various chicken lines is depicted in Table 2.
Overall, the numbers of upregulated genes are higher than downregulated genes, with line 15 having the highest number of unique DE genes (114 genes), and line O having the lowest unique DE genes (three genes).
DE genes shared by each of the chicken lines were illustrated in the representation of a six-tier Venn diagram (Fig. 3a, b) and the details of genes shared among the lines are shown in Table S6.The Venn diagram illustrated that 551 genes were commonly upregulated and 465 genes were commonly downregulated in all infected lines.In addition, 14 genes were identified to be upregulated in chicken lines with higher viral load (lines 6, 7, O and P; black circle, Fig. 3a), while 31 genes were upregulated in chicken lines with lower viral load (lines 15 and N; black circle, Fig. 3a).Several genes that were upregulated in lines 6, 7, O and P are CCR7, CST7 and VSIG4; while genes that were upregulated in lines 15 and N are DAPK2 and FAIM2 (Table S6).
On the other hand, two genes (NEDD8 and LOC428086) were identified to be downregulated in chicken lines 6, 7, O and P, while 44 genes were downregulated in chicken lines 15 and N (black circle, Fig. 3b).Genes that were downregulated in lines 15 and N are MAP2K6, IRF5 and TLR1.Genes that are associated with innate immune responses are TLR3, IL6, CXCLi2 and CCL4; though they were DE in the   infected chicken lines, their expressions are not differentially regulated by differences in the viral load.
The DE genes were further analysed by clustering genes corresponding to relevant functions, whereby each gene expression was compared between lines and viewed in a heat map (Fig. 4; detailed information of gene expression value is available in Table S7).The gene expression levels were clustered based on colour intensity, with red representing genes that were upregulated [log 2 (fold-change) >2)] and green representing genes that were downregulated [log 2 (foldchange) <À2] following vvIBDV infection.

Functional annotation and KEGG pathway enrichment analysis
The functional annotation and pathway enrichment analysis of 4588 DE genes were analysed based on the gene ontology (GO) and relevant KEGG pathways.A total of 28 KEGG pathways have been identified, where they were mainly involved in signalling pathways such as cytokine-cytokine receptor interaction, MAPK, JAK-STAT, Toll-like receptor (TLR), cell cycle, as well as replication and repair (Table S8).Nevertheless, the KEGG pathway enrichment analysis from the downregulated genes with log 2 (foldchange) <À2 reported that the genes were clustered into replication and repair, cell growth and death processes (Table S8).Meanwhile, a total of 57 GO terms were grouped into three ontologies; cellular component, molecular function and biological process (Fig. S2).

Expression of B-cell and T-cell-specific genes in the bursa following vvIBDV infection
BF which is a crucial organ in the B-cell development and differentiation showed downregulation of B-cell-specific genes including IGJ, CD79B, VPREB3, BLNK, TNFSF13B, TNFRSF13B, TNFRSF13C, EBF1, CD72 and IRF4; with the highest downregulation of VPREB3 gene in line 15 (Table S7), suggesting that all lines succumbed to IBDV infection.Nevertheless, upregulation of T-cell-specific genes, namely T-cell co-receptor CD3E, CD4, CD8 and LPAR (a G-protein coupled receptor induced in activated T-cells) were detected in most of the infected lines (Fig. 4).
Cytokines, chemokines and immune-specific genes responsive to vvIBDV infection Most of the genes involved in the immune responses were upregulated after the infection.Among them are IFNGR2, STAT1 and IRF1 that are involved in IFN-signalling pathways, and IFITM5, IFIT5, MX1 and SAMHD1, which are IFN-inducible genes.Expressions of several TLRs were also documented through the upregulation of TLR3 and TLR15 in all lines, downregulation of TLR2B in all lines and downregulation of TLR1A in lines 15 and N.Meanwhile, expressions of IFNB and IFNl3 (IL28B) were detected only in infected samples indicating that these genes were induced following vvIBDV infection.Other cytokine genes that were induced or upregulated in all lines were IL17A with the exception of lines 15 and O, and colony stimulating factor 1 (CSF1) except in line N.In addition, expression of CSF3 that was upregulated in lines N and P, was uniquely expressed in infected lines 15, 6, 7 and O (Fig. 4 and Table S7).
Several key pro-inflammatory cytokines and chemokines such as IL6, CXCLi2 (IL8) and CCL4 (MIP-1b) were upregulated in all lines.Meanwhile, genes such as IFNG, IL1B and CCL5 (RANTES) were upregulated in several lines.Although IL6 was highly expressed in line N at 9.96-fold, IL6ST (a critical component for IL6 signal transduction by forming a complex with IL6-specific receptor) was upregulated in all lines with the exception of line N.Meanwhile, CCR5, the receptor for CCL4 and CCL5, was found to be downregulated in line 15.Chemokine CXCL12, which was often induced by pro-inflammatory stimuli, was upregulated in line 7 at 2.91-fold; while its binding receptor CXCR4 was downregulated in both lines 6 and N. Other chemokines such as CCL19, CCL20, CCLI8, CCLI10, CCL14, XCL1 and CXCL13L2 were upregulated in most of the lines, while CXCR5, the binding receptor for CXCL13L2, was downregulated in all lines (Fig. 4).
IL12 is another pro-inflammatory cytokine whose expression can be detected based on its subunit IL12A or IL12B expression.IL12A was not DE in all lines; meanwhile expression of IL12B was differentially upregulated only in line N at 5.37-fold.In the meantime, IL15, which is closely related to IL2, was upregulated in both lines 7 and P; while its receptor IL15RA was upregulated in all lines except line 7.The DE genes that play an integral role in the host immune response as cytokines and chemokines are depicted in a cytokine-cytokine receptor interaction pathway (Fig. 5a).Other gene products implicated in the innate immune responses that were upregulated are as follows; GAL2, several proteases including lysozyme (LYG2), lysosomal cathepsin (CTSC, CTSD and CTSL2), matrix metalloproteinase (MMP2, MMP27, MMP28 and MMP23B) and granzymes (A and K).Moreover, upregulation of NOS2 was detected in all lines as expected due to macrophage activation.

Suppression of transcription and DNA replication factor following vvIBDV infection
Expressions of several genes involved in RNA transcription such as PAX3 and E2F8 were downregulated in all lines.Generally, the expressions of genes involved in DNA replication processes such as DNA2, FEN1, MCM-related genes, and RFC-related genes were also downregulated in every line (Fig. 4).

Stress response during vvIBDV infection
Alteration in genes associated with stress and unfolded protein response including those encoding for heat-shock protein HSPA2 (HSP70), HSPB1 (HSP27), HSPB7, HSPB8, HSP25, HSP90AB1 (HSP90) and SERPINH1 (HSP47) were detected in all lines.Expressions of HSPB1, HSPB8 and SERPINH1 were upregulated, while HSP90AB1 was downregulated.Meanwhile, the expression of HSPA2, which is often associated with virus infection, was upregulated in all   lines except line 7; with the highest upregulation in line O at 5.94-fold (Fig. 4).In addition to that, virus infection potentially stimulates expression of HSPA5 (GRP78), an unfolded protein response that is responsible for the activation of ER stress response [19].This gene was differentially upregulated only in line P (Fig. 4).

Apoptosis mechanism following vvIBDV infection
Several genes involved in the apoptosis cascades, such as pro-apoptotic genes CASP3 and DFFB, were downregulated in most of the lines with the highest downregulation of CASP3 observed in line 15 (Table S7).Anti-apoptotic gene, e.g.BCL2, which has been suggested to promote the survival of haemopoietic cells, was upregulated in lines 15 and 7. On the other hand, a member of the Bcl-2 family with proapoptotic function, BCL2L14 was upregulated in all lines with the exception of line N (Fig. 4).
In addition, the current transcriptomic profile of the expression of genes involved in p53 signalling showed a downregulated pattern including TP73 that regulates the transcription activity of p53, and also a wide spectrum of related genes such as CHEK2, GTSE, MDM4, STT3A, CDK1, CCNrelated genes, RRM2 and SESN3 (Fig. 4).This study also detected p38 MAPK activation through the upregulation of MAPK12 in all lines with the exception of line O, and MAPK13 in lines 15 and 7.However, the transcription factor MEF2C that is responsible for controlling cell proliferation and differentiation downstream of p38 MAPK was found to be downregulated in every line.Another important function of the p38 MAPK pathway is that its activation stimulates stress response by upregulating stress factor HSPB1 in all lines, with the highest upregulation exhibited in line 15 (Table S7).A pictorial illustration of the p38 MAPK signalling pathway is depicted in Fig. 5(b).

Gene expression study by real-time quantitative PCR
Expression of several genes with immune regulatory importance from the RNA-seq data was analysed using the realtime quantitative PCR platform, where relative expression of CASP3, CCL4, CXCLi2, IL12B, IFNG, TLR3, IRF10, BLNK and HSPA5 were expressed as log 2 fold-change.Meanwhile, expression of IL28B was presented based on corrected 40-cycle threshold (Ct) (Fig. 6a) since the expression of this gene is not detectable in mock-infected samples.Expression of IL28B on the real-time quantitative PCR platform was found to be the highest in lines 6 and 7 based on the corrected 40-Ct value (Fig. 6b).Overall, the real-time quantitative PCR approach detected slightly higher foldchanges than the RNA-seq platform for IL12B, IFNG, TLR3, CXCLi2, CCL4 and HSPA5.Nevertheless, the genes that were detected on the real-time quantitative PCR platform correlated with data generated by the RNA-seq platform with Pearson's correlation (R)>0.81,and the expression patterns of the selected genes were confirmed as previously detected by the RNA-seq platform (Fig. 6c).

DISCUSSION
Sequencing of the chicken genome has led to various studies in the description of numerous immune mediators including cytokines and chemokines that modulate immune responses towards infection [20,21].Inbred chickens with different major histocompatibility complex (MHC) haplotypes have been reported to demonstrate a different degree of susceptibility to disease towards invading pathogens [22].
In addition to the MHC gene, immune-related genes are also important in regulating host response against viral infection.Previous transcriptomic profiling studies were performed to determine host immune response following infection with classical IBDV strain [11,14].Therefore, there is a need to carry out transcriptional profiles of genes associated with immunological responses in different inbred chicken lines following infection with vvIBDV.
This study focuses on the use of the strand-specific NGS method to generate cDNA libraries for a robustly conserved transcript orientation and a more accurate determination of gene expression, as described in a previous study using mouse brain tissue [23].An average of 80 % annotated genes were successfully detected from read coverage of 50 M (79 bp) (Table S3) using the strand-specific sequencing approach.Additionally, the data generated in this study were validated using ERCC spike-in controls that served as external controls in all samples to ensure high-sequencing performance and an accurate fold-change threshold value measurement.A previous study on the bursal transcriptome highlighted that more genes were DE at 3 days p.i. compared to 1 days p.i. in SPF chickens infected with classical IBDV [16].In the current study, we explored genome-wide host responses of six different inbred lines based on the reference genome Gallus_gallus-4.0 to establish the bursal transcriptome baseline following vvIBDV infection, for further study on deciphering the roles of various immune mediators in modulating IBDV infection.Recently, a newly updated version of chicken genome is available at Gallus_gallus-5.0 (GCA_000002315.3).

It has been established that IBDV infection causes severe B-cell depletion and infiltration of T-cells together with
macrophages into the bursa of infected chickens [24,25].This result is in line with current findings, which stated that the decreased detection of transcription of B-cell-specific several were downregulated (orange).DE genes were also highly enriched in (b) the MAPK signalling pathway, indicating this pathway is crucial in modulating host immune response against vvIBDV infection.The green and white coloured boxes are default colours generated by the software, depicting genes that were not DE on RNA-seq, and unidentified genes in the organism-specific pathway, respectively.[4].In general, this finding suggests that vvIBDV infection caused moderate clinical signs [clinical score of (+5) >2 consecutive days], severe bursal lesion (score >3.5) and stimulated thousands of DE genes in all infected inbred lines, regardless of the differences in viral copy number at 3 days p.i.Further studies at a lower infectious dose (EID 50 ) of vvIBDV may be able to provide a more distinct differential expression of immune-mediated genes in the infected chickens.
The importance of the class I MHC gene during IBDV infection is also characterized since the inbred chicken lines used in this study expressed different class I MHC haplotypes [26].Line P with B19 haplotype, line O with B13 haplotype and lines 6 and 7 that shared B2 haplotype [26], showed higher viral copy number compared to line 15 with B15 haplotype and line N with B21 haplotype (Fig. 1).A previous study documented a link between the MHC haplotype with MHC class I expression through the BF2 gene, where expression of BF2 was lowest in line N, and highest in line P following MD infection across 40 chicken lines [27].In this study, upregulation of the BF2 gene was reported in lines P, 15, 6 and 7 (Fig. 4).Based on the viral copy number and bursal lesion scoring, line 6 of the B2 haplotype and line P of the B19 haplotype can be considered as more susceptible, while lines 15 and N as the least affected lines following vvIBDV infection.Unlike other lines, line O of the B13 haplotype, although showing the lowest bursal lesion score, had a high viral copy number, probably due to the overall lowest DE genes.Further study at different time points needs to be carried out to characterize the pathogenicity of this virus in these lines.
Infection with vvIBDV induced an upregulation of TLR3 and TLR15 (Fig. 4), where upregulation of TLR15 was also documented following MD virus and F52/70 IBDV infection [11,28].Other TLRs such as TLR2B and TLR1A were reported to be downregulated ensuing IBDV strain F52/70 infection [11].Another study has shown that chicken TLR1 and TLR2 are expressed in two isoforms A and B that formed heterodimeric complex interaction with different ligands, and influence host resistance towards pathogenic infection [29].In this study, we detected the downregulation of TLR2B in every line and the downregulation of TLR1A in lines 15 and N, which have lower viral load compared to other lines (Fig. 4).Therefore, future study is required to confirm the functional importance of TLR1 and TLR2 interaction during IBDV infection.
A new IFN-induced antiviral defence mechanism that has been documented in this study is type-III IFN, e.g.IFNl3 (IL28B).So far, no studies reported the IL28B activity in chicken, and we are the first to show IFNl3 expression following vvIBDV infection in vivo.IL28B was uniquely expressed since this gene can only be detected in the infected lines.A previous study has reported that expression of this gene can be induced by dsRNA in human foetal liver cells [30].The function of this gene in IBDV is unclear.However, a recent study showed that IFNl1 (IL29), a type-III chicken IFN, was able to inhibit in vitro influenza virus replication [31].
In this study, we reinforced the involvement of pro-inflammatory genes namely IL1B, IL6 and CXCLi2 in response to IBDV infection.The production of these pro-inflammatory cytokines and chemokines supports the involvement of macrophages, an important effector of innate immune response during IBDV infection [2].Upregulation of cytokines, chemokines and NOS2 indicates the occurrence of the inflammatory response in the infected bursa, which is in agreement with previous findings [32].Furthermore, studies have reported that an intense innate immune response especially the production of pro-inflammatory cytokines at an early phase of the infections may lead to severe bursal lesion and mortality, which may influence the host resistance mechanism depending on the virus strains and host's genetic background [7,11,14].In addition to that, the newly characterized avian CSF1 that plays an important role in regulating macrophage proliferation and differentiation [33], was upregulated in all the infected lines except line N (Fig. 4), indicating its involvement in macrophage activation following vvIBDV infection and subsequent promotion of the inflammatory response.
Chickens infected with IBDV generally showed an influx of CD4 + and CD8 + T-cell infiltration into the bursa, with the involvement of cytotoxic T-cells in the clearance of virusinfected cells through the granzyme-perforin pathway [34].
In this study, an upregulation of CD4 and CD8 expression was documented in all lines.Activation of T-cells in the bursa following vvIBDV infection is obvious through the expression of Th1 cytokines including IFNG and IL15, and also O-glycosylated transmembrane protein CD99.Another important Th1 cytokine involved in IBDV response was IL12B, an essential cytokine for IFNG expression suggesting vvIBDV infection in the bursa stimulates cell-mediated immune response, confirming earlier studies [35].Upregulation of granzyme A in lines 7, N, O and P, and NK-lysin in every line indicates the involvement of cytotoxic T-cells and/or NK cells in IBDV clearance, which was supported by earlier studies [11], stating that these cells were activated in the bursa following IBDV infection.
Expression of IFNs following virus infection activates the JAK-STAT signalling system that is important for the transcription of IFN-stimulated genes (ISGs), where ISGs are documented to be involved in antiviral defence by targeting any step in a virus life cycle to limit virus replication, and enhance IFN production [36].More than 100 avian ISGs have been successfully identified by RNA-seq and the microarray platform in the highly responsive chicken embryo fibroblasts (CEF) following type 1 IFN (IFNa) treatment [37].A previous study documented that recombinant chicken IFNa, and classical IBDV infection induced upregulation of several ISGs, including RSAD2, IFIT5, MX1 and EIF2AK2 (PKR), in CEF and in the bursa of SPF chickens, respectively [37,38].Similar findings were also documented in the bursa of vvIBDV-infected inbred lines, and other ISGs that were induced in this study include IFITM5, IRF7, IRF10, SAMHD1, SOCS3, TRIM25 and DHX58 (LGP2), suggesting the potentially important role of these ISGs in IBDV antiviral activities.One more important pathway involved in the pathogenesis of IBD is the p38 MAPK signalling pathway, i.e. it is an essential pathway to promote inflammatory response through macrophage activation that helps in the secretion of various cytokines and chemokines along with IL1B, IL6, CXCLi2 and NOS2, which corresponds with other studies [39].We also emphasize the involvement of the p38 MAPK pathway to induce a stress response following vvIBDV infection.It seems that the induction of pro-inflammatory cytokines in the MAPK pathway did not occur through the engagement of activator protein 1 (c-JUN) alone, but virus-induced inflammatory effect through secretion of IL1B.This later promotes the p38 MAPK pathway for the activation of HSPB1, a stress response gene required for downstream activation of proinflammatory mediators, with the ability to inhibit apoptosis by blocking CASP3 activity [40,41].Meanwhile, HSPA2 may also play a role in vvIBDV infection, where it has been reported to be upregulated in chicken brain tissue following H5N1 and Newcastle disease virus infection [42].
In this study, genes implicated in the apoptosis signalling pathway mostly were upregulated and involved as survival factors; including NGF, CSF2RB, BIRC2 and BCL2.Meanwhile, the genes that are involved in death factors were downregulated together with DFFB and CASP3.Downregulation of CASP3 was also recorded in a previous study of inbred chickens infected with F52/70 strain [14].Another pathway frequently implicated in IBDV-induced apoptosis is the p53 signalling pathway.A previous study showed that the p53 gene was upregulated at 2 days p.i., and downregulated later at 4 days p.i. in line 6 of IBDV-infected bursa [14].In the current study, a wide spectrum of related genes involved in p53 signalling were downregulated with earlier detection at 3 days p.i. following vvIBDV infection.Additionally, vvIBDV infection also stimulate the host immune responses through upregulation of genes involved in signalling pathways comprising MAPK, JAK-STAT, PPAR, TLR, TGF-beta and focal adhesion signalling pathways (Table S8).
In conclusion, all the different inbred lines showed positive immune response towards vvIBDV infection as they induced the activation of various mediators associated with macrophage and T-cell functions.Although some differences in the bursal lesion score and viral loads were detected in lines 7, 15, N and O, there was no clear explanation correlating this with the upregulation or downregulation pattern of gene expression.Further study using a lower dose of UK661 vvIBDV at different times p.i. at a much longer period (up to 7 days p.i.) may provide comprehensive variance in the clinical spectrum, as well as the susceptibility level of these inbred lines.The differential expressions of several mediators, namely IFNG, IL12B, IL15, IL28B, CCL5 and HSPB1 deserve further investigations as immunomodulator in the development of improved vaccine and therapeutic agent against IBD.S1.The birds were euthanized when the total clinical score reached (+5) for 2 consecutive days.By 3 days p.i., all the birds were euthanized.Bursa samples were collected from birds in each line of the control groups at 0 days p.i. (before virus inoculation) and from birds from each line of the infected groups at 3 days p.i.On the sampling day at 3 days p.i., the gross changes of the bursa were assessed and scored based on the scoring described by Shaw and Davidson [43].The remaining bursal tissues from each line was kept in RNAlater (Ambion, UK) solution and shipped to the Laboratory of Vaccines and Immunotherapeutics, Institute of Bioscience, UPM for NGS study.

RNA preparation
Bursal tissues that were collected from the infected groups (two birds per chicken line) and control groups (one bird per line) were used in the RNA-seq analysis.Each bursal tissue (20-30 mg) was homogenized using Tissue Ruptor (Qiagen, UK) in the presence of the RLT lysis buffer (Qiagen, UK).The total RNA was isolated using the RNeasy Plus Mini Kit (Qiagen, UK) according to the manufacturer's protocol.RNA concentration and purity was determined with the use of Bio-Spectrophotometer (Eppendorf, USA) and Bioanalyzer 2100 (Agilent Technologies, USA).Purified RNA samples at A260:A280 ratio !1.8 and RNA integrity number >8.0 were used for RNA sequencing.The prepared total RNA was spiked with 3-5 µl of 1 : 10 ERCC ExFold RNA Spike-in Control Mixes (Ambion, USA) serving as external RNA controls, where the sample mixtures were stored at À80 C until future use.

Preparation of cDNA libraries and read mapping
The purified RNA samples added with RNA spike-in were reverse-transcribed into cDNA libraries using TruSeq RNA library preparation kits, followed by standard Illumina pairedend cluster generation (Illumina, USA).DNA libraries with an insert size of ~200 nucleotides were clustered, amplified and sequenced on the Illumina HiSeq2000 platform according to the manufacturer's instruction.The sequence reads were assessed and pre-processed using FastQC [version 0.10.0 (www.bioinformatics.babraham.ac.uk/projects/fastqc/)] and FASTX-Toolkit [version 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/)] to obtain high-quality reads.The first 13 bases from the 5¢end, reads shorter than 30 bp with Q 20 or containing unknown bases were removed.The clean reads were mapped to the reference genome Gallus_gallus-4.0 (Gen-Bank accession ID: GCA_000002315.2), and also to the sequences of ERCC RNA spike-in, by a mapping tool TopHat (version 2.06) [44].

ERCC sensitivity assessment
The ERCC spike-in RNA controls were added to assess RNA-seq platform performances as described in the manufacturer's guideline.Briefly, the abundance of ERCC transcripts spiked in the samples was measured and normalized to fragments per kilobase of exon model per million mapped fragments (FPKM).The observed transcript's read density (log 2 FPKM) was plotted against its known transcript number, while the observed fold-change ratios for each ERCC control were plotted against the expected ratio to determine the dose response and fold-change response, respectively.

Transcript expression analysis
The alignment results of the reference genes were examined for transcript expression accompanied by differential expression analysis using Cufflinks and Cuffdiff [45], by allowing strand-specific processing.The expression level of the transcripts was expressed as FPKM, with a cut-off value of FPKM !1.0E-5 (determined from the ERCC doseresponse analysis).DE genes were identified with log 2 foldchange <À2 or log 2 fold-change >2 and Q-value <0.05.Genes were functionally clustered based on expression pattern using MeV: MultiExperiment Viewer [46].

Functional annotation and pathway enrichment analysis
The DE genes present in at least one chicken line were retrieved for further functional annotation and pathway enrichment analysis.KEGG pathways associated with the host response following vvIBDV infection were determined using DAVID Bioinformatic Resources [47], while functionally related DE genes were clustered with the corresponding GO terms using UniProt [48].

Real-time quantitative PCR assay
Total RNA samples used in Illumina were analysed with real-time quantitative PCR to validate the RNA-seq data.Briefly, 1 µg of RNA samples was reverse-transcribed into cDNA using the SensiFAST cDNA synthesis kit (Bioline, USA).Ten genes were selected with regard to their functions in the immune response and primer-probes were designed for the following genes: IFNG, IL12B, IL28B, TLR3, CCL4, CXCLi2, CASP3, BLNK, IRF10 and HSPA5, using TaqMan assays (Applied Biosystems, USA; Table S9).
Prior to quantification, total RNA extracted from bursal tissue infected with vvIBDV were used as a reference to generate a standard curve for each gene in a series of 10-fold dilution ranging from 1000 to 0.1 ng reaction on the CFX96 real-time system (Bio-Rad, USA).The data were expressed in Ct values to plot a standard curve to determine the optimized concentration of the template in the quantification experiment.For expression quantification, each cDNA sample was assayed in triplicate with the expression value normalized against two reference genes; ACTB and GAPDH, and the differential expression was expressed in fold-change from infected to mock-infected samples ratio using the -DDCT method [49].

Viral load quantification
IBDV gene-specific primers were designed as follows: forward primer; 5¢-ATG CTC CAG ATG GGG TAC TTC-3¢, and reverse primer; 5¢-TTG GAC CCG GTG TTC ACG-3¢ targeting the VP4 region of IBDV [50].The standard curve and melt curve were generated by SensiFAST SYBR kit (Bioline, USA) using the positive control cDNA of vvIBDV stock (10 7.4 EID 50 ml À1 ).A SYBR-based real-time quantitative PCR was performed in three-step cycling as follows: one cycle at 95 C for 2 min (polymerase activation), 40 cycles at 95 C for 5 s (denaturation), 65 C for 15 s (annealing step) and 72 C for 20 s (extension), while the final step of the melt curve analysis was performed at 70-95 C with increments of 0.5 C in every 5 s step À1 .All samples were assayed in triplicate and the Ct value was recorded to quantify the viral load based on log 10 transformation using the linear equation obtained from the standard curve generated at several different viral stock dilutions (10-fold dilution).

Statistical analysis
Expression of target genes and viral copy number quantification by real-time quantitative PCR were analysed for their significant value among lines using one-way ANOVA with Duncan's post hoc test where P<0.05 is considered as statistically significant.Expression data from RNA-seq was compared to expression data from real-time quantitative PCR using Pearson's correlation.

Fig. 1 .
Fig. 1.Viral copy number quantification.Viral copy numbers present in each sample were quantified based on absolute quantification by real-time quantitative PCR, and presented based on log 10 transformation of IBDV gene expression using the standard curve generated with slope: À3.569 and y-intercept: 60.045.The data shown are representative of two biological replicates and each biological replicate had three technical replicates.Statistical differences between inbred lines were assessed by one-way ANOVA (P<0.05).Means labelled with different letters are significantly different at P<0.05, and error bars represent the standard error of the mean.The figure bars showed that lines O, 7, 6 and P have significantly higher viral copy number compared to lines 15 and N.

Fig. 2 .
Fig. 2. Correlation analysis between replicates for inbred lines.Correlation analysis was performed for replicates of each inbred line, 15, 6, 7, N, O and P, and presented based on log 10 expressed transcripts (FPKM).The correlation coefficient (R 2 ) was between the range 0.83-0.95.

Fig. 3 .
Fig. 3. Identification of DE genes in different inbred lines at 3 days p.i. following vvIBDV infection.(a) Venn diagram depicting the uniquely expressed and differentially upregulated genes of the bursa collected from the vvIBDV-infected groups at 3 days p.i. compared to the mock-infected groups, with expression value log 2 (fold-change) >2 and Q-value <0.05.(b) Venn diagram depicting the genes that were absent and the differentially downregulated genes of the bursa collected from the vvIBDV-infected groups at 3 days p.i. compared to the mock-infected groups, with expression value log 2 (fold-change) <À2 and Q-value <0.05.The corresponding coloured-region represents each chicken line as follows; green: line 15, red: line 6, yellow: line 7, pink: line N, purple: line O and light blue: line P. Numbers displayed in each region indicate the numbers of common genes in different chicken lines.For example, 14 genes (circle) in (a) were common among lines 6, 7, O and P as upregulated genes.

Fig. 4 .
Fig. 4. Transcriptional response of chicken genes towards vvIBDV infection.Gene expression profiles of the bursa infected with UK661 vvIBDV at 3 days p.i. obtained from the RNA-seq was clustered according to functional categories for each inbred line (lines 15, 6, 7, N, O and P).The heat maps indicate the log 2 (fold-change) of the bursal gene expression with >2 (red) are differentially upregulated or <À2 (green) are differentially downregulated.Genes whose expression was not DE with log 2 (fold-change) between À2 and 2 are coloured in black.

Fig. 5 .
Fig.5.KEGG pathway enrichment analysis of the host response towards vvIBDV infection in the bursa.These KEGG pathways were generated using DE genes from all inbred lines (lines 15, 6, 7, N, O and P) that were obtained from RNA-seq.Most of the DE genes were highly enriched in (a) the cytokine-cytokine receptor interaction pathway with most of them being upregulated (blue) whilst, only

Farhanah 35 Fig. 6 .
Fig. 6.RNA-seq platform validation by real-time quantitative PCR.(a) Relative expression of genes selected from RNA-seq data, namely CASP3, CCL4, CXCLi2, IL12B, IFNG, TLR3, IRF10, BLNK, HSPA5 and IL28B were measured on the real-time quantitative PCR platform based on the -DDCT method.The vertical axis represents log 2 fold-changes of gene expression.The vertical axis for IL28 represents corrected 40-Ct value since this gene was absent from the mock-infected samples.Data represented as mean log 2 foldchange ±SEM and mean of corrected 40-Ct±SEM of triplicate measurements.Statistical difference between inbred lines were assessed by one-way ANOVA (P<0.05).Means labelled with different letters are significantly different at P<0.05, and error bars represent the SEM.(b) Changes in IL28B expression following vvIBDV infection was compared between RNA-seq and real-time quantitative PCR platforms.IL28B expression was reported based on the FPKM value of infected samples using the RNA-seq platform, meanwhile, the expression on the real-time quantitative PCR platform was expressed as corrected 40-Ct.(c) Correlation of the gene expression level between real-time quantitative PCR and RNA-seq platforms.The expression level for nine genes (IFNG, IL12B, TLR3, CCL4, CXCLi2, CASP3, BLNK, IRF10 and HSPA5) from the bursal-infected vvIBDV were analysed for their correlation based on log 2 fold-change between the two platforms for each chicken line (lines 15, 6, 7, N, O and P).The dots represent the genes and the R in the graphs indicates the correlation coefficient.** represents the significant level [P<0.01;using Pearson's correlation for linear relationship (twotailed)].Downloaded from www.microbiologyresearch.orgby IP: 54.70.40.11On: Sun, 30 Dec 2018 08:09:22

Table 1 .
Clinical signs and bursal lesion score in six different inbred lines following vvIBDV infection

Table 2 .
Summary of DE genes in inbred chicken lines following vvIBDV infectionAll genes mapped to the reference genome on RNA-seq were analysed for differential expression with threshold values log 2 (fold-change) >2 or log 2 (fold-change) <À2 and Q-value <0.05.
*DE genes only in the infected samples considered as uniquely expressed.†Absent in infected samples/expression in infected samples or did not pass the detection threshold level.Downloaded from www.microbiologyresearch.orgby IP: 54.70.40.11On: Sun, 30 Dec 2018 08:09:22 be associated with the decreased number of Bcells.Based on the viral copy number and bursal lesion score, line 6 has the highest bursal lesion score and a high viral copy number, while line 15 has a low bursal lesion score and a low viral copy number (Fig.1and Table1).In contrast, line O has the lowest average bursal lesion score and a high viral copy number.Nevertheless, it is interesting to note that line O has the lowest DE genes.The susceptibil- [12]loaded from www.microbiologyresearch.orgbyIP:54.70.40.11On:Sun, 30 Dec 2018 08:09:22genes might ity level of vvIBDV-infected inbred lines was previously determined by Bumstead et al.[12], based on the mortality rate that generally occurred between 3-8 days p.i.[12].Due to ethical constrains in the current study, further changes in clinical development and mortality could not be carried out beyond 3 days p.i., and therefore susceptibility level of these inbred lines could not be determined exclusively.The use of a single time point (3 days p.i.) to profile transcriptome expression across six inbred lines provides an unequivocal insight into the host response during the clinical endpoint of the disease, since previous study has shown severe clinical disease caused by UK661 vvIBDV observed between 56 and 72 h p.i. in RIR chickens Inbred SPF chicken lines; 15, 6, 7, N, O and P from unvaccinated flocks were supplied by the Institute of Animal Health, Compton.The percentage heterozygosities of the 7) and 13.3 % (15).These inbred lines have been maintained through over 20 generations by full-sibling mating.The inbred lines with mixed gender were randomly assigned into two experimental groups; infected and control groups.Infected and age-matched control birds were housed in separate experimental animal rooms and supplied with both water and a vegetable-based feed ad libitum.