Experimental challenge with bovine respiratory syncytial virus in dairy calves: bronchial lymph node transcriptome response

Bovine Respiratory Disease (BRD) is the leading cause of mortality in calves. The objective of this study was to examine the response of the host’s bronchial lymph node transcriptome to Bovine Respiratory Syncytial Virus (BRSV) in a controlled viral challenge. Holstein-Friesian calves were either inoculated with virus (103.5 TCID50/ml × 15 ml) (n = 12) or mock challenged with phosphate buffered saline (n = 6). Clinical signs were scored daily and blood was collected for haematology counts, until euthanasia at day 7 post-challenge. RNA was extracted and sequenced (75 bp paired-end) from bronchial lymph nodes. Sequence reads were aligned to the UMD3.1 bovine reference genome and differential gene expression analysis was performed using EdgeR. There was a clear separation between BRSV challenged and control calves based on gene expression changes, despite an observed mild clinical manifestation of the disease. Therefore, measuring host gene expression levels may be beneficial for the diagnosis of subclinical BRD. There were 934 differentially expressed genes (DEG) (p < 0.05, FDR <0.1, fold change >2) between the BRSV challenged and control calves. Over-represented gene ontology terms, pathways and molecular functions, among the DEG, were associated with immune responses. The top enriched pathways included interferon signaling, granzyme B signaling and pathogen pattern recognition receptors, which are responsible for the cytotoxic responses necessary to eliminate the virus.

Bovine Respiratory Disease (BRD) is the leading cause of morbidity and mortality in calves over one month of age in Ireland 1-4 and internationally [5][6][7][8] . It is a multifactorial disease, involving infectious agents (both viral and bacterial), host factors, management practices, environmental stress factors and their interactions 4,[8][9][10][11] . Viruses generally initiate BRD and the damage they cause to the respiratory tract predisposes the calves to a secondary infection caused by bacteria, many of which are normal flora of the bovine upper respiratory tract 12 .
Bovine Respiratory Syncytial Virus (BRSV) is one of the main viral infectious agents responsible for the onset of BRD 13 . BRSV is an enveloped, non-segmented, negative-stranded RNA virus and is a member of the Orthopneumovirus genus within the family Pneumoviridae 14,15 . The morbidity resulting from infection with BRSV is in the range of 60% to 80%, and mortality is reported to be 20% 15 . The disease resulting from infection with BRSV can vary from sub-clinical to severe manifestation of clinical signs including coughing, pyrexia, nasal discharge, anorexia and respiratory distress 15 . Lung pathology resulting from BRSV infection is primarily due to the host inflammatory response, particularly the induction of pro-inflammatory cytokines (IL12B, IFNG, TNF, IL6, IL18, CXCL8, CCL5, CCL2, CCL3, IFNA1 and IFNB1) and the influx of leukocytes, particularly neutrophils 13,15 . Establishment and maintenance of BRSV infection is facilitated through the virus' ability to interfere with the host's interferon (anti-viral) response and to induce immunomodulation by shifting the adaptive immune response away from a cell mediated response towards a Th2 dominated response 15 www.nature.com/scientificreports www.nature.com/scientificreports/ Animal sampling. On days −1, 0, 6 and 7, relative to the challenge (day 0), a 9 ml K 3 EDTA blood tube was collected, gently inverted several times and placed on ice (Fig. 1). Whole blood was analysed for haematological variables (white blood cell count, neutrophil percentage, lymphocyte percentage, monocyte percentage and eosinophil percentage) on an ABAXIS H5 haematology analyser (ABRAXIS Model: H5 S/N: 364372) immediately following collection. Daily, from the challenge to the day of slaughter, clinical signs (nasal discharge, ocular discharge, general demeanour, size of mandibular lymph nodes, presence of a cough, respiratory rate, respiratory character, mouth breathing, dyspnoea, presence of an expiratory grunt and rectal temperature) were recorded and scored by a veterinarian, who was blinded to the calves' treatment status (BRSV challenged or control), (Fig. 1) using a modified version of the clinical scoring system described by Gershwin et al. 21 (Supplementary Table S1). Using this scoring system, points were accrued for each abnormal clinical symptom and the total number of points corresponded with the severity of disease, such that higher clinical scores were associated with more severe BRD 21 .
On the day of the challenge, calves were either administered with 10 3.5 TCID 50 /ml × 15 ml of BRSV inoculum (BRSV challenged calves in groups B and C) or mock challenged with 10 ml PBS (control calves in group A) by aerosol inhalation.
Animals were restrained in a cattle crush and the head was further restrained to prevent movement during nebulisation. Animals were nebulised with BRSV strain SVA 274/92 20 or PBS using the Omron NE-U780 nebuliser (OMRON Model: NE-U780-E S/N: 20151200182AF) fitted with a veterinary face mask (GaleMed Model: VM-2 Animal Mask size 5 22mmlD Ref: 5635 Lot: 151019). Airflow and nebulisation volume settings were set to the maximum. BRSV inoculum (20 ml) was placed into the chamber and nebulised over 15 minutes. The remaining volume of inoculum was checked after 15 minutes and if less than 10 ml had been nebulised, nebulisation was continued for an additional 5-10 minutes.
On day 7 relative to the challenge, calves were euthanized by captive bolt. The lungs were scored for lesions by a veterinarian using an AFBI standard lung scoring system which assigns the percentages of lesions present on the total lung area and on component parts of the lungs ( Supplementary Fig. S1).
The workspace and implements were thoroughly cleaned and disinfected with bleach and 75% ethanol and sprayed with RNaseZap for RNase inhibition, before tissue collection and between euthanisation of animals. Bronchial lymph node tissues were harvested immediately and flash frozen in liquid nitrogen, placed on dry ice and then transferred to a −80 °C freezer for storage.

Statistical analysis of clinical and haematological data.
Clinical score and haematological variables (white blood cell number, neutrophil percentage, lymphocyte percentage, monocyte percentage and eosinophil percentage) were analysed using repeated measures mixed models (MIXED procedure of SAS v 9.4) where time-point defined the repeated measure. Data were first assessed for normality using PROC REG and PROC UNIVARIATE. All data were found to be normally distributed. Treatment (BRSV challenged or control), time-point or group (A, B or C) and their interactions were included as fixed effects. Calf was included as a random effect. A Tukey adjustment was used to correct for multiple testing. Means and their standard errors are presented in Figs 2 and S2 (graphs were generated in GraphPad Prism version 7.02). Lung scores were assessed for normality using PROC REG and PROC UNIVARIATE and analysed using a mixed model ANOVA (MIXED procedure of SAS v 9.4) with treatment (BRSV challenged or control) or group included as fixed effects. The presence or absence of lung lesions was analysed using a Fisher's Exact test in SAS 9.4.

RNA extraction.
Total RNA (including miRNAs) was extracted using the Qiagen RNeasy Plus Universal Mini Kit (Qiagen LTD, Manchester, UK), according to the manufacturer's instructions including Appendix C of the manufacturer's protocol. The quantity of the extracted RNA was determined by measuring the absorbance at 260 nm with a Nanodrop spectrophotometer (NanoDrop technologies; Wilmington, DE, USA). The quality of the extracted RNA was examined with the Agilent 2100 Bioanalyser (Agilent Technologies Ireland Ltd; Dublin, Ireland) using the RNA 6000 Nano LabChip kit (Agilent Technologies Ireland Ltd; Dublin, Ireland). Samples had a mean RNA Integrity Number of 8.6 ± 0.31. www.nature.com/scientificreports www.nature.com/scientificreports/ RNA-Seq library preparation and sequencing. Extracted RNA was shipped frozen at −80 °C on dry ice to the University of Missouri's DNA Core Facility for RNA-Seq library preparation using the TruSeq stranded mRNA Kit (Illumina, San Diego, California, USA) and high-throughput sequencing (75 bp paired-end) on an Illumina NextSeq 500. All sequence data produced in this study has been deposited to NCBI GEO repository and are available through the series accession number GSE131452.
Alignment of sequence reads to the bovine reference genome and differential gene expression analysis. Adapter trimmed sequence reads in FASTQ format were assessed for quality using FastQC (version 0.11.7) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). All reads passed the basic quality statistics. Reads were aligned to the UMD3.1 bovine reference genome and read counts were generated by converting aligned reads into counts per gene using the Spliced Transcripts Alignment to a Reference (STAR) aligner (version 2.5.2b) 22 .
Differential gene expression was determined using the R (R version 3.4.2 2017-09-28) Bioconductor package EdgeR (version 3.20.9) 23 which accounts for biological and technical variation by modelling data as a negative binomial distribution. To remove lowly expressed genes, any genes with less than one count per million in at least three of the samples, were removed from the analysis. Data were normalised across libraries using the trimmed mean of M-values normalisation method 24 . The quantile-adjusted conditional maximum likelihood (qCML) common dispersion and the qCML tagwise dispersions were used to estimate dispersion. Exact tests were used for the detection of DEG between BRSV challenged and control calves. Genes with a Benjamini-Hochberg false discovery rate (FDR) of 10% and a fold change of ≥2 were considered differentially expressed. pathway and gene ontology analysis. Statistically significant gene ontology terms, pathways and other gene function related terms were analysed from a ranked DEG list (DEG between BRSV challenged and control calves with an FDR of 10% and a fold change of >2, ranked by P-value) using g:Profiler version 2018-06-21 25 . Multiple testing correction was performed using the g:Profiler tailor-made algorithm g:SCS.
The DEG between BRSV challenged and control calves, with an FDR of 10% and a fold change of >2, were analysed for over-represented gene ontology terms/pathways and functional annotation clustering using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/tools.jsp) 26 . The enrichment P-values for each annotation term were derived from a modified Fisher's exact test called the EASE score. The group enrichment score for each functional annotation cluster was calculated from the geometric mean (in the -log scale) of all the enrichment P-values in each annotation term within each cluster and was used to determine the significance of results from the functional annotation clustering analysis. A group enrichment score threshold of 1.3 (EASE = 0.05) was applied.
To further examine over-represented pathways, cellular and molecular functions and predicted upstream regulators, the RNA-Seq data were analysed using the Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., https://www. qiagenbioinformatics.com/products/ingenuitypathway-analysis), according to the manufacturer's instructions 27 . Within IPA, Fisher's exact test was used with the Benjamini-Hochberg correction for multiple testing for the identification of over-represented pathways and over-represented molecular and cellular functions with a FDR of 10%, from DEG between BRSV challenged and control calves, with an FDR of 10% and a fold change of >2. Additionally, IPA's regulation Z-score algorithm, which predicts increases or decreases in functions based on directional changes in the DEG and expectations derived from the literature, was used to predict differences in the over-represented cellular and molecular functions. IPA software considered cellular and molecular functions with a regulation Z-score value of ≥2.0 to be significantly increased and cellular and molecular functions with a regulation Z-score value of ≤−2.0 to be significantly decreased.

Comparison to the BRSV-associated RNA-Seq Data from Tizioto et al. (2015). FASTQ files were
downloaded from the National Center for Biotechnology Information Sequence Read Archive under accession number SRP052314 (SRR1952354, SRR1952355, SRR1952356, SRR1952357, SRR1952358, SRR1952359, SRR1952360) 18 . Alignment, gene counts, differential gene expression analysis, pathway and gene ontology analyses were performed identically to the analyses described for the current study. Pathway and gene ontology analyses were additionally performed on genes which were commonly differentially expressed in this study and that of Tizioto et al. 18 with an FDR of 10% and a fold change of >2.

Results
clinical score. Clinical score was affected by treatment (BRSV challenged 42.0 ± = 4.02 versus control 27.3 ± 5.68; P = 0.03) (Fig. 2). Time-point of clinical assessment also affected clinical scores (P < 0.0001) ( Fig. 2) with days 4 (52.7 ± 7.42) and 7 (45.2 ± 7.42) being significantly greater than day 0 (25.5 ± 7.42). There was a tendency towards an interaction between time-point and treatment affecting clinical scores (P = 0.1), with the BRSV challenged calves having numerically greater clinical scores on day 7 versus day 2 (P < 0.1) and the control calves having significantly greater clinical scores on days 4 and 6 compared to day 2 (P < 0.05). However, there were no significant differences or tendencies between the BRSV challenged calves and the control calves at any of the timepoints (P > 0.05). Furthermore, there was no effect of group on clinical scores (P > 0.05).
Lung pathology. Figure 3 shows the lung pathologies from a subset of the BRSV challenged calves. There were no differences in either the overall lung scores or the percentage of the right cranial lung lobe that was lesioned between the BRSV challenged calves and the controls (P > 0.05). This is likely due to the small percentages of lesioned tissues present on the lungs (BRSV challenged group average percentage of lesioned tissue overall was 1.34 ± 1.0% and the right cranial lobe 1.63 ± 0.8%, controls 0.08 ± 0.1% and 0.83 ± 0.8%) ( Table 1).
However, overall, only one control calf exhibited a lesioned lung, while 5 out of the 6 calves in BRSV challenge group B and 4 out of the 6 calves in BRSV challenge group C had lung lesions ( Table 1). The BRSV challenged calves had a higher probability of having a lesioned lung compared with the control calves (P = 0.04).
Sequence alignment. An average of 53,878,694 (2 × 75 bp) paired-end reads (i.e. 26,939,347 sequenced fragments) were generated for each sample (Supplementary Table S2). Overall, 7.2% of the sequence reads failed to map to the genome and 89.6% of sequence reads mapped uniquely to the genome. On average, 57.3% of the total number of sequence reads were assigned to a gene by STAR, 0.4% were ambiguous and 31.9% of the sequence reads aligned to no feature (Supplementary Table S2).
Differential gene expression and functional annotation. Using Edge-R, a multi-dimensional scaling (MDS) plot was produced which allowed an examination of the similarity of the samples based on global gene expression when samples were plotted using their first two principal components. The MDS plot indicated no differences between individuals in BRSV challenge groups B and C, but a clear separation between the BRSV challenged and the control calves (Fig. 4).  Table S3). There were two DEG between the calves in BRSV challenge group B and BRSV challenge group C (PPP1R1B and LYNX1).
From the g:Profiler analysis, there were 133 over-represented biological process Gene Ontology terms, 26 over-represented cellular component Gene Ontology terms, 23 over-represented molecular function Gene Ontology terms, 17 molecular function KEGG pathways and 20 over-represented Reactome pathways (P corrected <0.05) (Supplementary Table S4). The 12 most significant over-represented biological process Gene Ontology terms were related to the immune response ( Table 2). The over-represented KEGG pathways were also associated with immune response and diseases and included KEGG:05164 Influenza A (Fig. 5) 28 . This pathway contained 19 up-regulated and two down-regulated genes (Supplementary Table S4, Fig. 5).
There were 77 over-represented ontology terms in the DAVID analysis (p < 0.05, Benjamini-Hochberg FDR <0.1) (Supplementary Table S5). These were associated with the immune response, glycoproteins and secreted peptides. There were 28 enriched functional annotation clusters in the DAVID analysis (group enrichment score >1.3) (Supplementary Table S6). The enriched clusters contained terms relating to signal peptides, immunity,    Table 2. The 12 most significant over-represented biological process Gene Ontology terms among the genes which were differentially expressed between BRSV challenged and control calves. www.nature.com/scientificreports www.nature.com/scientificreports/ immunoglobulin, nucleotide-binding, thrombospondin, endoplasmic reticulum, double-stranded RNA binding, viral diseases, inflammatory responses, lectins, chemokine activity, extracellular matrix structure, cell division, calcium binding and the complement pathway.
Ingenuity pathway analysis showed that seventeen canonical pathways were enriched (P < 0.05, FDR <0.1) (Supplementary Table S7, Fig. 6). Four of these pathways were predicted to be upregulated (Interferon Signaling, Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses, Inhibition of Matrix Metalloproteases and Acute Phase Response Signaling) (Z-score >2) (Supplementary Table S7, Fig. 6).
Thirty-one enriched biological functions were predicted to be increased in BRSV challenged animals compared to controls (Z-score >2) (Supplementary Table S8). Their functions were related to the recruitment, activation and death of leukocytes.
There were 256 up-stream regulators that were predicted by IPA to be activated (Z-score ≥2) (Supplementary Table S9). Nine of these had an observed increase in fold change (IFNG, IRF7, DDX58, CGAS, OSM, IL27, EIF2AK2, XBP1 and TBX21). There were 84 up-stream regulators that were predicted to be inhibited by IPA (Z-score ≤2) (Supplementary Table S9), but none of these were found to have a decreased expression in the BRSV challenged calves. were obtained from each BRSV challenged/control bronchial lymph node sample downloaded from the National Center for Biotechnology Information Sequence Read Archive under accession number SRP052314 (Supplementary  Table S2) 18 . Overall, 11.7% of sequence reads failed to map to the genome and 82.7% of reads mapped uniquely to the genome, using the STAR aligner 22 . On average, 28.0% of the total number of sequence reads were assigned to a gene by STAR, 0.1% were ambiguous and 54.6% of reads aligned to no feature (Supplementary Table S2).

Differential gene expression and functional annotation of RNA-Seq data in comparison to
An MDS plot illustrating the similarity of the Tizioto et al. 18 BRSV challenged and control bronchial lymph node samples based on log 2 fold gene expression changes was produced in EdgeR (Fig. 7). For illustration of the similarities between the samples in the present study and those from the Tizioto et al. 18 study, an MDS plot was produced in EdgeR based on log 2 fold gene expression changes among all samples (Fig. 8).
One thousand and thirteen genes from the study of Tizioto et al. 18 were differentially expressed between BRSV challenged and control calves (Supplementary Table S10). From the g:Profiler analysis, there were 198 over-represented biological process Gene Ontology terms, 19 over-represented cellular component Gene Ontology terms, 16 over-represented molecular function Gene Ontology terms, 18 over-represented KEGG pathways, 38 over-represented Human Phenotype Ontology terms, 1 over-represented CORUM term and 30 over-represented Reactome pathways (P corrected <0.05) (Supplementary Table S11). The over-represented biological process Gene Ontology terms were related to the immune response. The over-represented KEGG pathways were also associated with immune responses and diseases and included KEGG:05164 Influenza A and KEGG:04060 Cytokine-cytokine receptor interaction.
There were 354 DEG in common between the present study and the Tizioto et al. 18 study (BRSV challenged and control bronchial lymph node samples) (Supplementary Table S12). From the g:Profiler analysis, there were 81 over-represented biological process Gene Ontology terms, 18 over-represented cellular component Gene Ontology terms, 13 over-represented molecular function Gene Ontology terms, 16 over-represented KEGG pathways, 12 over-represented Human Phenotype Ontology terms, 1 over-represented CORUM term and 20 over-represented Reactome pathways (P corrected <0.05) among the DEG in common between the two data sets (Supplementary Table S13). www.nature.com/scientificreports www.nature.com/scientificreports/ Pathway analysis in IPA revealed that 15 pathways were over-represented among the DEG in common between the present study and the Tizioto et al. 18 study (Fig. 9, Supplementary Table S14). The most significant of the over-represented pathways, interferon signaling and the role of pattern recognition receptors in recognition of bacteria and viruses, were also among the most significant over-represented pathways in the IPA pathway analysis of the DEG solely found in this study.
Fourteen enriched biological functions, among DEG in common between the present study and the Tizioto et al. 18 study, were predicted to be increased in BRSV challenged animals compared to controls (Z-score >2) (Supplementary Table S15). Functions were related to the antimicrobial and antiviral immune response. Twelve enriched biological functions, among DEG in common between the present study and the Tizioto et al. 18 study, were predicted to be decreased in BRSV challenged animals compared to controls (Z-score >2) (Supplementary  Table S15). These functions were associated with viral replication and morbidity.

Discussion
This is the first study to examine the host response at the transcriptional level in artificially-reared Holstein-Friesian calves to a challenge with BRSV. Previous studies have examined the bronchial lymph node pathological response to a BRSV challenge in Angus-Hereford crossbred beef calves 18,19,21 . However, immune responses, particularly at the level of transcription, have been demonstrated to be affected by breed [29][30][31] . This  www.nature.com/scientificreports www.nature.com/scientificreports/ may be the reason that the majority of the variation observed in the MDS plot (Fig. 8) was due to study (the present study versus Tizioto et al. 18 ), while treatment (BRSV challenged versus control) bronchial lymph node gene expression was the secondary source of variation. However, the different sequencing platforms used (NextSeq 500 versus HiSeq 2000), the varying length and depth of reads generated, the distinct BRSV viral strains used in the inoculum and the alternative husbandry and management practices carried out in the studies' diverse geographical locations, may have also influenced the alterations in gene expression between the two studies. Despite the many sources of variation existing between the two studies, 38% of the DEG in this study were commonly differentially expressed in the RNA-Seq analysis of the BRSV challenged and control bronchial lymph node samples from the study of Tizioto et al. 18 . Furthermore, the majority of enriched pathways found in the two studies were in common. Therefore, BRSV appears to elicit many conserved alterations in gene expression.
In the present study, calves were challenged with a similar titre of BRSV that has consistently and successfully been used to establish infection in previous BRSV challenge studies 18,21,32 . Calves were euthanised 7 days post challenge as this is the time-point when the acute phase response has been observed to peak 33 and when live BRSV and its antigen titres begin to decrease 34 . Bronchial lymph node tissues were harvested and RNA-Seq analysis was performed on these tissues to determine the gene expression changes induced by BRSV infection as these lymph nodes drain the lungs and are the site of antigen presentation and the initiation of the primary immune response 35,36 .
Infection with BRSV can be subclinical, present only in the upper respiratory tract or it can affect both the upper and lower parts of the respiratory tract 15 . BRSV mainly infects and produces clinical signs in calves less than one year of age but is usually more severe in young calves aged between one and three months 32 . It is common for paramyxoviruses, including BRSV, to be attenuated by passage through in vitro culture, and this can result in minimal clinical symptoms of respiratory disease 37 . Therefore, the age of the calves (approximately 5 months of age) and the effects of BRSV attenuation in culture may have influenced the primarily sub-clinical response to BRSV challenge observed in the present study. The calves in this study had substantially lower clinical scores than the calves challenged with BRSV in the Tizioto et al. 18,21 study. This may have been due to the possible differing initial BRSV specific MDA levels of the calves, the diverse breeds of the calves or the distinct strains of BRSV inoculum used for the challenge. Furthermore, there were no changes in neutrophil, lymphocyte or monocyte percentages due to the BRSV challenge in this study. While control calves did have a higher white blood cell count compared with BRSV challenged calves throughout the study, the white blood cell numbers were within the normal range for healthy cattle 38 and no increase in white blood cell count was observed following virus inoculation, in BRSV challenged calves.
Interestingly, despite these BRSV challenged calves only showing a modest clinical and haematological response to infection, and having only a small proportion of their lungs lesioned, the virus induced vast changes in bronchial lymph node gene expression levels. Therefore, measuring host gene expression levels may prove to be a more accurate way of diagnosing BRD than by the observation of clinical signs of the disease. This would be hugely beneficial to the global agricultural industry as sub-clinical BRD which remains undiagnosed is almost as prevalent in healthy cattle as BRD is among cattle which are diagnosed with respiratory disease and treated for BRD 39 . Although it would clearly not be possible to measure host gene expression in bronchial lymph node tissue in live calves if a small set of diagnostic genes are also perturbed in expression in white blood cells or nasal fluid, these more accessible media could be used as "liquid biopsies" for the future diagnosis of subclinical respiratory infections. Indeed increased expression of IFNG and IL6 have been observed along with clinical signs during an www.nature.com/scientificreports www.nature.com/scientificreports/ initial BRSV challenge and the increased expression of IFNG in the peripheral blood mononuclear cells of young calves with an absence of clinical symptoms has also been observed in a reinfection with BRSV 34 .
The majority of the alterations in bronchial lymph node gene expression due to BRSV challenge detected in the present study were associated with the immune response. The most highly upregulated gene was GZMB, which was upregulated sixty-seven fold. This gene was also the most up-regulated in response to a BRSV challenge in bronchial lymph node in US calves 18 and the encoded protein has been shown to be significantly increased in CD8+ and CD4+ T cells and natural killer cells, in tracheal aspirates from children infected with respiratory syncytial virus (RSV) 40 . Furthermore, GZMB was found to be a member of several enriched ontologies and pathways including "defense response", "immune response", "granzyme-mediated apoptotic signaling pathway" and "granzyme B signaling", in the present study. The serine protease protein encoded by this gene is one of the major components of the lytic granules which are released by cytotoxic T cells and natural killer cells to destroy virally infected cells 41 . It also plays a role in degradation of the extracellular matrix and induction of pro-inflammatory cytokine responses 42 . Furthermore, granzyme B has recently been found to be utilised by T regulatory cells to control effector B and T cell populations, thereby regulating the inflammatory response and reducing the severity of RSV disease in mice 43 . Consequently, it is likely to play an important role in defence to BRSV infection and alterations to this gene may be critically important to BRD resistance and susceptibility. Indeed, it has been shown that polymorphisms in murine GZMB can reduce the ability of CD8+ T cells to kill virally infected target cells 44 . It has also been postulated that that GZMB-expressing T regulatory cells could be used as a therapeutic strategy in human viral lung infections, due to their ability to dampen inflammation without elevating viral load 43 . However, a study of mice experimentally infected with Klebsiella pneumoniae found that while GZMB knockout mice had higher levels of lung inflammation at earlier stages in the disease, there were minimal differences in distant organ injury or survival between the knockout and wild type mice 45 . However, Klebsiella pneumoniae is caused by a bacterium rather than a virus and as GZMB is mainly involved in the cell mediated cytotoxic immune response, it may be more effective in the elimination of viral infections and intracellular bacteria.
Pathways and gene ontologies involved in the innate immune response, induction of cytokines and chemokines, humoral and cell mediated immunity were over-represented among the DEG. The KEGG pathway for Influenza A was enriched for DEG and this pathway contains upregulated cytokines responsible for the induction of the febrile response, mitochondrial genes involved in programmed cell death and interferon stimulated genes which induce cytotoxic T cell responses. Tizioto et al. 18 also observed an enrichment of similar pathways and the DEG found in common in both studies were members of the most significantly enriched pathways from this study, including interferon signaling and role of pattern recognition receptors in recognition of bacteria and viruses. Furthermore, the induction of pro-inflammatory chemokines and cytokines in calves due to BRSV infection has been reported in pneumonic lung lesions 46 , bronchoalveolar lung fluid 47 and bovine respiratory/turbinate epithelial cells 48 . Therefore, these immune response associated genes seem to play an important role in the bovine host response to BRSV infection. This is expected as pathogen recognition, induction of inflammatory cytokines and chemokines, and the initiation of the innate immune response are vital to control the virus upon entry. Subsequent elimination of the virus requires the development of CD8 cytotoxic T cells, which is facilitated by the production of IFNG and interferon stimulated genes such as DDX58, IFI44, IFIH1, IFIT1, IFIT2, IFIT3, IFIT5, ISG15, ISG20, MX1, MX2, OAS1,  EIF2AK2, RTP4, RSAD2, OAS2 and OSM, which were upregulated in the present study 49 .
In agreement with Tizioto et al. 18 , we also predicted neurodegenerative diseases and functions including "progressive motor neuropathy", "progressive neurological disorder", "inflammatory demyelinating disease" and "multiple sclerosis" to be increased in the BRSV challenged calves. These findings are interesting as human RSV has been observed as being capable of infecting murine lung neuronal cells 50 , and also of infecting the human central nervous system and inducing neurological symptoms 51 .
BRSV modulates the immune response towards a T helper 2 (Th2) type response, with increased expression of Th2 promoting cytokines and increased concentrations of BRSV specific IgE antibodies in lymph fluid 48 . Despite enriched ontology terms including "B cell activation" and "immunoglobulin mediated immune response", we found no evidence for an increase in the transcription of either IL4 or IL13 which initiate the Th2 bias 52 . This is consistent with the results of Grell et al. 34 , who did not detect an increase in IL4 gene expression following a BRSV challenge. However, although there was a 2 fold increase in the expression of IL12RB2, a transmembrane protein subunit of the interleukin 12 receptor complex, there was no increase in the expression of IL12, which is necessary to induce Th1 cell proliferation responses 48 , which is consistent with the expected Th2 bias for BRSV. Lipopolysaccharide (LPS) was a predicted upstream regulator in this study and that of Tizioto et al. 18 . Furthermore, LBP was the second greatest up-regulated gene in response to the BRSV challenge infection in this study (fold change = 52). This is likely because the fusion protein from BRSV interacts with TLR4 to initiate the innate immune response in a manner that is similar to LPS 48 . Tizioto et al. 18 suggested that LPS has the potential to be used in new treatments aiming to prevent the occurrence of BRD as it induces a cascade of both anti-viral and anti-bacterial genes. Other predicted upstream regulators in the present study include IFNG and IRF7, which play roles in IFN antiviral innate immune responses 53 .
BRSV is closely related to human RSV, and induces similar pathologies 15,37 and alterations in host gene expression. Genes involved in the type 1 interferon antiviral response, including BP1, RSAD2, MX2, MX1, ISG15, IRF7, IFITM3, IFIT3, IFIT2, IFIT1, IFIH1, IFI44, which were upregulated in the lungs of C57BL/6 mice, in response to infection with RSV 54 , were also up-regulated in the calves challenged with BRSV in the present study. Therefore, the gene expression changes induced by BRSV infection detected in this study may provide an insight into the human bronchial lymph node transcriptional response to RSV. Future studies will also elucidate the bovine host micro-RNA gene expression responses and determine the regulatory regions controlling these transcriptional alterations due to infection with BRSV, in bronchial lymph node, and in other important respiratory tissues including lesioned and healthy lung tissue, mediastinal and retropharyngeal lymph nodes, and the pharyngeal tonsil. The DEG identified in this study, particularly those also observed to be DEG in the study of Tizioto et al. 18 , and therefore involved in