Transcriptomic analysis among sub-clinically ill cattle following Bovine Respiratory Disease vaccine protocol and challenge with Bovine Viral Diarrhea Virus

Abstract

Background

Bovine viral diarrhea (BVD) is a widespread disease of cattle affecting cattle production on both a domestic and global basis [1, 2]. Bovine viral diarrhea virus (BVDV) causes an infection that can often go unnoticed and yet be persistent in herds. Health impacts on infected cattle range from mild to mortal, and include diarrhea, respiratory complications, mucosal diseases, abortions, and hemorrhagic complications [3-5].

Complicating the issue with BVDV infections is the virus predisposes both infected and non-infected cattle to bovine respiratory disease complex (BRD) that can lead to widespread contagions in herds [6]. Losses in the cattle industry due to BVDV infected herds and BRD in general have increased costs to maintain production and treat animals, and reduced outputs of beef or dairy programs. In spite of the focus and effort to understand and increase the immune response to both BVDV and BRD pathogens, there are gaps in the knowledge of vaccine efficacy, whether there are differences between the protection afforded by modified live (MLV) or killed (KV) vaccines used to safeguard herds against BRD, and if they offer cross-protection to BVDV. 

There also exists a gap in the knowledge of the genomic elements involved in the efficacy of a particular vaccine strategy when signs of infection are sub-clinical. The issue with sub-clinical infections can be exacerbated when feedlot interactions are coupled with the high pathogenicity of BVDV infections, controlling the spread of the virus as a preventative measure becomes a daunting task. This prompts the need for research to explore the genetic underpinnings of tolerance that when coupled with proper herd vaccinations, which could lead to a more robust method of controlling the spread of BVDV even when infections are undetectable because they are sub-clinical. Research into this topic has the additional benefit of providing producers with information on vaccine efficacy that can help them adopt vaccination strategies that reduce treatment costs.

The objective of the study was to determine the extent of cross-protection by examining the underlying differential gene expression response to a sub-clinical BVDV challenge in MLV and KV BRD vaccinated Nellore-Angus cattle.

Results

Clinical data

Analysis of the antibody titers at 14 dpi revealed that the MLV group had lower titers (P < 0.01) for IBR, BVD1a, BVD1b, and BVD2 when compared to the KV group (Fig. 1). There were no detectable antibody titers in the non-vaccinated (NON) group at 14 days post infection (dpi).

Transcriptomic analysis shows difference in immune arm activation based on vaccine type

Clear immunosuppressive effects of the BVDV challenge were observed in all group comparisons. Both the MLV and KV groups showed downregulation of key host anti-viral genes such as MX dynamin like GTPase 2 (MX2) and 2'-5'-oligoadenylate synthetase 2 (OAS2) [7, 8]. However, the analysis also showed that immunosuppression was either directed towards the adaptive immune response as seen with the KV group or the innate immune response as seen in the MLV response. What is unclear from the study is the effect of time on these responses as only one time point (14 dpi) was observed. The lateness of this sample date could explain a portion of the muted immune responses observed within the study, and it is likely that the majority of the responses are linked to the ability of BVDV to subdue the host immune response as late as 14 dpi [4, 5].  With this in mind, our results indicated that the BRD MLV vaccine offered better cross-protection to ensuing BVDV infection supported by underlying differential expression of multiple immune genes involved in MHC, T-cell, intestinal, and adaptive immunity. However, both vaccine strategies revealed possible candidate genes that could help elucidate the differences in cross-protection experienced by the animals (Table 1). Results of the transcriptomic analysis allowed us to examine two phenomena within the study: the first being the mRNA response to the difference in cross-protection by vaccine type. The second being the immunosuppressed transcriptomic response that underlies sub-clinical BVDV infections.

Differential expression between killed vaccine and non-vaccinated groups

Comparison of expression between the BRD KV and NON groups resulted in 38 differentially expressed genes (FDR < 0.15) (Table S1). Most of these genes were found to be downregulated (n= 29) and skewed towards adaptive immune functions [7, 8]. This included downregulation of immune genes such as joining chain of multimeric IgA and IgM (JCHAIN) involved in the biological processes of the adaptive immune response and leukocyte migration. This also included the observed downregulation of anti-viral genes such as MX2, OAS2, and interleukin 23 receptor (IL23R). Of the few upregulated genes that could be identified, multidrug resistance-associated protein 4-like (LOC509854) had the highest fold-change (log2fc = 1.25).  This gene and stood out because of functions related to the KEGG pathways for bile secretion and anti-folate resistance. The anti-folate resistance pathway is an anti-inflammatory pathway found to be active during the administration of drugs meant to control microbial infections and chronic inflammation brought on by disease in livestock. At first glance the connection to the antifolate resistance pathways appears to extraneous, however research has been conducted into a group of antifolate compounds named Artemisinins that were shown to limit replication of BVDV virions, therefore behaving anti-virally [9, 10].  Other upregulated genes of interest included ATP binding cassette subfamily C member 4 (ABCC4), another multi-drug resistance gene also involved in cellular redox processes, and LOC618416 leukocyte immunoglobulin-like receptor subfamily A member 6, an inhibitory immune function gene [11], one of the only adaptive immune response genes shown to be overexpressed in the KV group. Downstream protein-protein interaction prediction and G.O. analysis of the KV and NON group further revealed that inhibition of adaptive immune function pathways was present at 14 dpi. This inhibition is a possible hallmark of the KV vaccines ability to be less protective observable through statistically significant (FDR < 0.1) G.O. results.

Differential expression between modified live vaccine and non-vaccinated groups

Comparison of the  MLV and NON groups revealed 82 statistically significant DE genes, which although more than between the KV and NON group, showed the same propensity to be skewed toward downregulation (16 upregulated/66 downregulated) (Table S2). Some of the downregulated genes in this comparison have been shown to exert positive effects on the adaptive arm of the immunity such as the gene Immunoglobulin superfamily member 2 (CD101;IGSF2), which is an empirical inhibitor of T-cell proliferation [7, 8, 12]. It may be that downregulation of CD101 facilitates the added protection of the MLV vaccine. 

Several upregulated immune genes of interest with adaptive immune actions appeared in the MLV vs. NON contrast. This included the CD1d molecule (CD1d) that facilitates lipid and glycolipid antigen presentation to T cells. Additionally, CD1d participates in the KEGG pathways for viral process, the innate immune response, and the adaptive immune response [7, 8, 13]. The other upregulated genes of interest were enteric beta defensin (EBD; DEFB1); a member of a group of neutrophil derived microbicidal and cytotoxic bovine beta defensins [7, 8, 14] and the major histocompatibility complex, class I, A (BoLa). It is possible that the combined overexpression of these several genes is supporting the lower MLV titers and improving MLV cross-protection. The downstream gene ontology analysis revealed downregulation of multiple anti-viral reactome pathways that included ISG15 antiviral mechanism and antiviral mechanism by IFN-stimulated genes (FDR = 1.49E-02; pathway log2fc = -0.73). This result is likely driven by the overwhelming number of downregulated immune genes appearing within the list of differentially expressed genes for the MLV vs. NON contrast.

Differential expression between modified live and killed vaccine groups

The final contrast was between the two vaccinated sample groups and contained the largest observable number of differentially expressed genes (n = 186) (Table S3), wherein 41 of the genes were upregulated, and 145 were downregulated. The large number of statistically significant genes in this list spotlights the differences in efficacy observed between the MLV and KV groups. Overall there was mostly downregulated immune genes in the MLV group when contrasted against the KV group with most of the genes known to operate as part of the bovine anti-viral immunity. This implies that the KV group had mostly upregulated innate immune genes despite having higher antibody titers. The up- and downregulated genes of interest overlapped with some those from the MLV vs. NON contrast, such as EBD, IL23R, and CD101, however the magnitude of expression observed was greater. Additionally, more upregulated genes related to triggering innate and cell-adaptive immunity were observed for the MLV group. This included the major histocompatibility complex, class II, DQ alpha 5 (BOLA-DQA5) involved in the intestinal immune network for IgA production KEGG pathway; complement C1q C chain (C1QC) and complement C1q B chain (C1QB) both part of the initial triggering of classical complement cascade portion of the innate immune response [7, 8].

The differences in the MLV and KV cross-protective effect were the most evident during the examination of the gene ontology groupings based on the sub-clinical transcriptomic response (Table 2). Statistically significant (FDR < 0.1) immune processes for mononuclear cell proliferation (GO:0032943; FDR = 3.71E-02); Th1 and Th2 cell differentiation (KEGG:04658; FDR = 2.63E-02) were composed of mainly downregulated genes, while the gene ontology groupings for T-cell activation (GO:0042110; FDR = 1.07E-02) was composed of both up- and downregulated genes. The Reactome analysis revealed more downregulated immune pathways within the MLV vs. KV contrast that appear to show a dulling of harmful inflammatory responses such as neutrophil degranulation (FDR = 1.38E-04; pathway log2fc =-0.82), an immune response that when prolonged can cause tissue damage, but nonetheless is a marker of reduced innate immune responses in the MLV group [15]. Other downregulated reactome pathways included the innate immune system (FDR = 3.67E-03; pathway log2fc = -0.65) and the ECM proteoglycans (FDR = 3.67E-03; pathway log2fc =-0.96), which could be a response to the ability of certain BVDV strains to lead to mucosal disease development [16]. There was only one upregulated Reactome pathway observed, the classical antibody-mediated complement activation pathway (FDR = 9.17E-02; pathway log2fc = 0.83), which is involved in the innate immune response.  

Discussion

The immune response to pathogens is an amalgamation of the host’s innate and adaptive immunity. The innate immune response exists as a broad range of non-specific primary responses to infection, while the adaptive immune response represents acquired immunity. Part of the work of the innate immune response allows for the priming of adaptive immunity within a host. A branch of the innate immune response, thought to have a role in bridging the innate and adaptive immune responses [17, 18],  is represented by the complement activation pathways which acts to recognize pathogens and then help eliminate them through proteolytic chain of reactions that lead to proinflammatory stimulation, lysis, or opsonization that “complements” the activity of the innate immune cells increasing the host’s defensive response to pathogens. There are several complement pathways that bolster the innate immune response through sequential events that lead to activation of C3 and C5 convertase; the classical pathway which is initiated by antigen-antibody recognition, the lectin pathway initiated by pathogen recognition by lectins, and the alternative pathway [19]. Within this study, the classical complement pathway was shown to be active in the MLV vaccinated cattle. This pathway bolsters the innate and adaptive immune response and is initiated by the recognition of pathogen surfaces and /or antibody-antigen complexes formed as part of the adaptive immune response by the gene C1q. The Gene C1q is the first step in the classical complement pathway cascade. Structurally, functional C1q is composed of several genes; C1QA, C1QB, and C1QC, which contain globular domains that allow for binding to the antigen-antibody complexes. Initiation of the classical pathway through C1q leads to a bolstering of the immune response through pro-inflammatory peptides and phagocytosis of pathogens [18, 20]. The components of C1q; C1QA, C1QB and C1QC were observed to be overexpressed (FDR ≤ 0.15) within the MLV group in this study, but not the KV group. This shows a potential difference in the how each vaccine reinforces innate and adaptive immune functions. The ability of the classical complement pathway to recognize the antigen-antibody complexes allows it to act as a bridge between the innate and adaptive responses that furthers the effectiveness of the combined immune response.

The results of the transcriptomic analysis of the sub-clinical BVDV infection uncovered multiple genes and pathways related to differences in the immune response and cross-protection related to BRD vaccination. This study allowed for inroads into the immune expression that may accompany sub-clinical BVDV infections, however larger sample sizes will be needed to fully explore this phenomenon. The challenge presented by BVDV and BRD infections portends a need for identifying ways for cattle producers to find the best level of protection balanced against the cost of treatment. While the KV may be considered cheaper and safer, the MLV has the advantage of not needing a booster and appears to lead to immune system differences at the genomic level.

Increased complement activation in MLV group may bolster better cross-protective immune response

It is possible that within the MLV vaccinated cattle, the classical pathway of the complement system is being stimulated by the over-expression of C1QA, C1QB and C1QC, which make up functional C1q, to maintain primed phagocytic cells that can bridge the innate and adaptive immune response to the viral cells. This overexpression may allow for quicker pathogen recognition and immunoregulatory activation [18] through the complement cascade not initiated in the KV vaccinated cattle, in which C1QC was downregulated. It could also be feasible that the upregulation of the complement cascade system is a side effect of the immunosuppressive activity exerted by BVDV during infections and may be related to negative regulation by C1QC and C1QB of host macrophages and granulocytes [7, 8, 20, 21]. Cai et al. [22] showed that significantly upregulated C1q can be a potential biomarker of active tuberculosis infection [22]. It may be that the slight upregulation in the samples marks a slow return to the baseline expression levels after adaptive immunity initiation at 14 dpi, which may be useful as a marker of sub-clinical infection. The MLV group may also be experiencing a better response through the added complement pathway mediated pathogen recognition and neutralization.

Overexpression of key genes related to promoting adaptive immunity were unique to the MLV vaccinated group

The immune related genes were overexpressed in the MLV treated group included the complement cascade genes; C1QA, C1QB, and C1QC, discussed previously and other immune genes such as MHC gene BOLA-DQA5, EBD, CD1d, and a cytokine sequestration gene in IL23R.

The gene IL23R may be the key to the MLV group having a better response to the BVDV challenge due to being necessary to stimulate CD4⁺ maturation into Th17 pro-inflammatory molecules [23].  Because no statistically significant differential expression of the interleukin 23 gene was observed, it is possible that IL23R overexpression at 14 dpi is a remnant of the biological processes involved in Th17 production that helped drive the MLV group toward adaptive immunity by creating a more robust innate immune response. Additionally, IL23R is a key player in intestinal immune responses and has been dubbed an inflammatory bowel disease (IBD) susceptibility gene [24].

The gene BOLA-DQA5 is one of the major bovine MHC-II genes [25] and indicated that the MLV group had transitioned into adaptive immune functions despite lower antibody titers. Besides antigen presentation, the gene also takes part in Th1, 2, and 17 differentiation which may help bolster MLV efficacy in concert with some of the other genes observed within the study. The gene is also a component of the intestinal immune network for IgA production and its over expression could help limit mucosal disease prevalence in MLV vaccinated cattle.

The gene EBD is a bovine beta defensin gene involved in IL-17 signaling. In general, the beta-defensins function as antimicrobial peptides and chemo-attractants in cattle [26, 27]. However, EBD is considered an enteric beta defensin involved in intestinal immunity and may serve as an indication of the importance of long-term neutrophil potentiation in the protection against BVDV. The overexpression of this gene in intestinal tissue of cattle could possibly serve as a marker of adaptive protection unique to MLV vaccination.

The gene CD1d functions as an antigen presenting glycoprotein that helps drive innate immune functions against viral pathogens through presentation to natural killer T-cells and the MHC. The gene is considered to be MHC-like and involved in the immunoregulatory action of lymphoid cells. It is possible that the overexpression of CD1d is linked to higher mean levels of lymphocytes observed in the MLV group at 14 dpi, which may have allowed for better viral clearance by the innate immune system effectively reducing the volume of neutralizing antibodies needed by adaptive immune system. Another possibility is that CD1d overexpression led to increased antigen presentation that effectively increased the adaptive response to react earlier and led to less antibody production at later timepoints. Interestingly, studies in mice have shown viruses such as the herpes simplex virus (HSV) can inhibit CD1d antigen presentation [13] that could assist viral evasion of the host immune system. However, it seems that in the MLV treated group CD1d has the ability to escape the immunosuppressive effects of BVDV.

The difference in immune class potentiation between MLV (adaptive) and KV (innate) could explain the difference in anamnestic antibody response

Downey-Slinker et al. [28] showed evidence that the increased antibody response observed within the KV group compared to the MLV was not a reflection of protection, but rather a reflection of lower MLV titres, which may be more indicative of reduced viral replication within the MLV treated group. Downey-Slinker et al. [28] also suggested that it is possible that this difference is linked to more than just a humoral response. This was observable within our study as the MLV group displayed a higher number of cell-mediated immune response genes being differentially expressed in response to challenge. Although the KV group had higher antibody titres than the MLV group, it did not have many upregulated adaptive immune response genes being expressed when contrasted against the MLV group. Also observable in the MLV vs. KV contrast was that the KV group showed more innate immune gene upregulation at 14 dpi, suggesting the KV group was still actively battling the virus. In the KV treatment group one of the most upregulated genes is the LOC618416 leukocyte immunoglobulin-like receptor subfamily A member 6 gene, which is the bovine homolog for the human LILRB5 gene. These receptors are found on the surface of antigen-presenting cells. The human homolog of LOC618416 is considered to have immunoregulatory functions and may elucidate the lack of adaptive immune gene expression in the KV group. In humans these receptors can bind to MHC genes through immunoreceptor tyrosine-based inhibitory motifs that effectively inhibit signaling of adaptive immune responses [7, 11, 29].

It is possible that upregulation of LOC618416 is preventing the antigen-presenting cells in the KV group from properly signaling the MHC, leading to little or no adaptive immunity being triggered during sub-clinical BVDV infections. This could be the result of the immunosuppressive ability of BVDV. Although, the KV group has been shown to have higher antibody titers at 14 dpi for all challenge pathogens, the differential expression analysis showed that many of the innate immune function genes were still being upregulated. This provides additional evidence that the KV group had not fully transitioned away from the innate response by 14 dpi in contrast to the MLV group.

Interestingly, some of the human inhibitory leukocyte immunoglobulin-like receptors have been shown to interact with the gene CD1d [30], which was only seen as being upregulated in the MLV vs. KV and MLV vs. NON comparisons. It may be possible that part of the difference in vaccine efficacy is linked to upregulated leukocyte immunoglobulin-like receptors impairing the antigen presenting activity of CD1d that could possibly delay onset of adaptive immunity [13].

What may also contribute to the difference between MLV and KV cross-protective efficacy is the downregulation of CD101 observed in the MLV-treated group. In humans, expression of CD101 inhibits normal T-cell proliferation effecting both lymphocyte activation and IL2 pro-inflammatory signaling [12, 31]. The study provides evidence that part of the difference in vaccine efficacy of the MLV vs. KV treated groups can be represented by the differential expression of genes involved in lymphatic activity and antigen presentation to T-cells. The upregulation of CD1d, EBD, and IL23R along with the downregulation of CD101 within the MLV group appears to allow for a more efficient activation of the adaptive immunity possibly through better T-cell immunomodulation. This was observable through the expressed MHC and MHC-like genes at 14 dpi that were only seen to be upregulated in the MLV contrasts. This activity at 14 dpi within the MLV group is juxtaposed to the upregulation of genes with immune inhibition functions witnessed in the KV group.

Conclusions

Both vaccines provide some level of cross-protection, however the results point to MLV producing cross-protection that extends the adaptive immune response to BVDV. Overall, the upregulation of MHC molecules and T-cell related gene expression observed within the MLV comparisons indicate that the BRD-MLV supports cross-protective immune expression that dulls the immunosuppressive effect that some BVDV strains exert on adaptive immunity. The current study provides researchers with useful information on the connection between vaccine efficacy and host immune responses, but is only a first step. Future research can expand on this study’s results by examining additional BVDV strains and increasing the sample size to better evaluate the immune response in sub-clinical infections.

Methods

Animals and sample collection

Animal selection, vaccination, and challenge procedures are explained in full in Runyan et al.  [32]. In brief, the cattle used in the study were part of the Texas A and M University McGregor Genomics herd and the animal protocol was approved by the Texas A&M University Institutional Animal Care and Use Committee and the Texas A&M University Institutional Biosafety Committee. The project examines gene expression of white blood cells from three treatment groups; no vaccine (NON), killed vaccine (KV), and modified-Live vaccine (MLV) of F₃ Nellore-Angus crossbred cattle intranasally infected with a type 1b non-cytopathic BVDV [2] after vaccination with commercially available BRD vaccines. The killed vaccine (KV) groups received an initial vaccination of Vira Shield^® (Novartis Animal Health US, Inc., Greensboro, NC, USA) at -56 days and a booster at -35 days prior to challenge [32].  The MLV group received a single vaccination at day -35 of Arsenal^®4.1 (Novartis Animal Health US, Inc., Greensboro, NC, USA), and were kept isolated from KV and non-vaccinated animals for 7 days [32]. At day 0 animals were administered 5 ml (1 × 10⁵ TCID/mL) of BVDV1b strain CA0401186a intranasally [2]. Leukocytes were separated as buffy coats within 60 min of sample collection 14 dpi and snap frozen in liquid nitrogen. Additional phenotypic observations of clinical signs of infection were recorded each day from 1 to 14 dpi. Only information from 14 dpi was used to evaluate and filter samples for animals showing only sub-clinical signs of infection for downstream analysis.

RNA isolation and sequencing library preparation

Libraries for RNA sequencing were generated from RNA extracted from bovine leukocytes using the TruSeq Total RNA Library RiboZero prep kit ^® (Illumina, Inc., San Diego, CA). Two chips on an Illumina NextSeq 500 were used to generate 150 bp single-end reads for 18 libraries.

Transcriptomic analysis

Sample sequences were concatenated and checked for quality using FASTQC prior to trimming based on quality score (> 32) and read length that discarded any reads shorter than 18 bps. Quality control measures resulted in a range of ~50-60 million reads per sample being retained for alignment to the UMD 3.1 bovine reference genome using the HiSAT software package within the Galaxy bioinformatics portal [33, 34]. The hybrid nature of the Nellore (Bos indicus)/Angus (Bos taurus) samples contributed to an average total (unique and multiple) sample alignment of 80-86%. The samples were then filtered based on collected clinical phenotype data indicating which samples were negative of visible clinical signs. This filter excluded 1 sample from each treatment group (MLV, KV, NON) leaving a total of five samples per vaccine treatment group to be analyzed for differential gene expression using DESeq2. The model used to examine sub-clinical gene expression of the leukocytes included variables for ~ vaccine treatment + rectal temperature + average daily weight gain at 14 dpi + lungscore + pen. Annotation of gene counts was performed using the FeatureCounts software package [35]. Differential gene expression was calculated using the DeSeq2 package with the dispersion model set as local and all other parameters set at their default values [36]. Statistical significance was set to FDR ≤0.15 to account for the nature of sub-clinical responses, which are hypothesized to result in more subtle changes in transcript expression.

Pathway and gene ontology analysis

Downstream analysis of gene biological pathways was carried out using a combination of online gene ontology and over-enrichment software packages that included Reactome, STRING 10.5, and g:GOST [37-39]. Reactome analysis used Ensembl gene IDs for Bos taurus and did not project to human or use other interactors; analysis using the g:GOST software was based on the default parameters except the size of query/term intersection parameter was changed to reflect pathways containing at least 3 query terms from the input list; and the STRING 10.5 analysis was based on the parameters of molecular action, query terms only, and a medium confidence interaction score (> 0.400). The pathway and gene ontology analysis used a statistical significance threshold of FDR ≤ 0.15 based on the Benjamini and Hochberg FDR adjustment for multiple gene set corrections [40].

Statistical analysis

Statistical analysis of the comparison of anamnestic antibody response for each treatment group was conducted using ANOVA followed by Tukey-Kramer HSD test in JMP^® (Version 13 SAS Institute Inc., Cary, NC, 1989-2007). The anamnestic antibody response levels for each group were calculated by the method described in Downey-Slinker et al. [28].

Abbreviations

BVDV:  Bovine viral diarrhea virus

BVD:  Bovine viral diarrhea

BRD:  Bovine respiratory disease

MLV:  Modified-live vaccine

KV:  Killed vaccine

NON:  Non-vaccinated

DPI: Days post infection

Declarations

Ethics approval and consent to participate

The animal use protocol was reviewed and approved by the Institutional Animal Care and

Use Committee (IACUC) of Texas A and M University.

Consent for publication

Not applicable

Availability of data and material

All relevant data are within the paper and its Supporting Information files. Data files

have been submitted to the GEO database ( )

Competing interests

CAG is a member of the editorial board (Associate Editor) of this journal. The remaining authors declare they have no competing interests.

Funding

This work was supported by the NIFA AFRI grant 2017-67012-25999 from the USDA National Institute of Food and Agriculture. Funding for collection of animal data was provided by Texas A&M AgriLife Research through its Beef Competitiveness Program and was a component of USDA CSREES Hatch project TEX08937. The funding body had no role in the design of the study and collection, analysis and interpretation of data, or in writing the manuscript.

Authors’ contributions

DSF and CAG designed the transcriptomics study.  ADH led the vaccination and challenge study and coordinated phenotype collection.  CAG extracted the leukocytes.  DSF coordinated RNAseq, performed all data analyses, and drafted the first version of the manuscript under the supervision of CAG.  LCM assisted with approach to analysis and manuscript revisions. All the authors read and approved the final manuscript.

Acknowledgements

We would like to thank the members of the Gill and Herring lab groups for the technical assistance in collecting samples and underlying data. We would also like to thank Dr. Andrew Hillhouse and TIGGS for technical assistance in sequence preparation.

References

1. Ridpath JF: Bovine viral diarrhea virus: global status. Vet Clin North Am Food Anim Pract 2010, 26(1):105-121, table of contents.
2. Ridpath JF, Neill JD, Peterhans E: Impact of variation in acute virulence of BVDV1 strains on design of better vaccine efficacy challenge models. Vaccine 2007, 25(47):8058-8066.
3. Ridpath J: The contribution of infections with bovine viral diarrhea viruses to bovine respiratory disease. Vet Clin North Am Food Anim Pract 2010, 26(2):335-348.
4. Chase CC: The impact of BVDV infection on adaptive immunity. Biologicals 2013, 41(1):52-60.
5. Peterhans E, Schweizer M: BVDV: a pestivirus inducing tolerance of the innate immune response. Biologicals 2013, 41(1):39-51.
6. Larson RL: Bovine Viral Diarrhea Virus-Associated Disease in Feedlot Cattle. Vet Clin North Am Food Anim Pract 2015, 31(3):367-380, vi.
7. Jenuth JP: The NCBI. Publicly available tools and resources on the Web. Methods Mol Biol 2000, 132:301-312.
8. The UniProt C: UniProt: the universal protein knowledgebase. Nucleic Acids Res 2017, 45(D1):D158-D169.
9. Romero MR, Serrano MA, Vallejo M, Efferth T, Alvarez M, Marin JJ: Antiviral effect of artemisinin from Artemisia annua against a model member of the Flaviviridae family, the bovine viral diarrhoea virus (BVDV). Planta Med 2006, 72(13):1169-1174.
10. Krishna S, Bustamante L, Haynes RK, Staines HM: Artemisinins: their growing importance in medicine. Trends Pharmacol Sci 2008, 29(10):520-527.
11. Hogan L, Bhuju S, Jones DC, Laing K, Trowsdale J, Butcher P, Singh M, Vordermeier M, Allen RL: Characterisation of bovine leukocyte Ig-like receptors. PLoS One 2012, 7(4):e34291.
12. Soares LR, Tsavaler L, Rivas A, Engleman EG: V7 (CD101) ligation inhibits TCR/CD3-induced IL-2 production by blocking Ca2+ flux and nuclear factor of activated T cell nuclear translocation. J Immunol 1998, 161(1):209-217.
13. Horst D, Geerdink RJ, Gram AM, Stoppelenburg AJ, Ressing ME: Hiding lipid presentation: viral interference with CD1d-restricted invariant natural killer T (iNKT) cell activation. Viruses 2012, 4(10):2379-2399.
14. Meade KG, Cormican P, Narciandi F, Lloyd A, O'Farrelly C: Bovine beta-defensin gene family: opportunities to improve animal health? Physiol Genomics 2014, 46(1):17-28.
15. Lacy P: Mechanisms of degranulation in neutrophils. Allergy Asthma Clin Immunol 2006, 2(3):98-108.
16. Shimshoni E, Yablecovitch D, Baram L, Dotan I, Sagi I: ECM remodelling in IBD: innocent bystander or partner in crime? The emerging role of extracellular molecular events in sustaining intestinal inflammation. Gut 2015, 64(3):367-372.
17. Dunkelberger JR, Song WC: Complement and its role in innate and adaptive immune responses. Cell Res 2010, 20(1):34-50.
18. Sontheimer RD, Racila E, Racila DM: C1q: its functions within the innate and adaptive immune responses and its role in lupus autoimmunity. J Invest Dermatol 2005, 125(1):14-23.
19. Sarma JV, Ward PA: The complement system. Cell Tissue Res 2011, 343(1):227-235.
20. Lu JH, Teh BK, Wang L, Wang YN, Tan YS, Lai MC, Reid KB: The classical and regulatory functions of C1q in immunity and autoimmunity. Cell Mol Immunol 2008, 5(1):9-21.
21. Nayak A, Pednekar L, Reid KB, Kishore U: Complement and non-complement activating functions of C1q: a prototypical innate immune molecule. Innate Immun 2012, 18(2):350-363.
22. Cai Y, Yang Q, Tang Y, Zhang M, Liu H, Zhang G, Deng Q, Huang J, Gao Z, Zhou B et al: Increased complement C1q level marks active disease in human tuberculosis. PLoS One 2014, 9(3):e92340.
23. Cauli A, Piga M, Floris A, Mathieu A: Current perspective on the role of the interleukin-23/interleukin-17 axis in inflammation and disease (chronic arthritis and psoriasis). Immunotargets Ther 2015, 4:185-190.
24. Pidasheva S, Trifari S, Phillips A, Hackney JA, Ma Y, Smith A, Sohn SJ, Spits H, Little RD, Behrens TW et al: Functional studies on the IBD susceptibility gene IL23R implicate reduced receptor function in the protective genetic variant R381Q. PLoS One 2011, 6(10):e25038.
25. Gelhaus A, Forster B, Wippern C, Horstmann RD: Evidence for an additional cattle DQA locus, BoLA-DQA5. Immunogenetics 1999, 49(4):321-327.
26. Roosen S, Exner K, Paul S, Schroder JM, Kalm E, Looft C: Bovine beta-defensins: identification and characterization of novel bovine beta-defensin genes and their expression in mammary gland tissue. Mamm Genome 2004, 15(10):834-842.
27. Tarver AP, Clark DP, Diamond G, Russell JP, Erdjument-Bromage H, Tempst P, Cohen KS, Jones DE, Sweeney RW, Wines M et al: Enteric beta-defensin: molecular cloning and characterization of a gene with inducible intestinal epithelial cell expression associated with Cryptosporidium parvum infection. Infect Immun 1998, 66(3):1045-1056.
28. Downey-Slinker ED, Ridpath JF, Sawyer JE, Skow LC, Herring AD: Antibody titers to vaccination are not predictive of level of protection against a BVDV type 1b challenge in Bos indicus - Bos taurus steers. Vaccine 2016, 34(42):5053-5059.
29. Brown D, Trowsdale J, Allen R: The LILR family: modulators of innate and adaptive immune pathways in health and disease. Tissue Antigens 2004, 64(3):215-225.
30. Burshtyn DN, Morcos C: The Expanding Spectrum of Ligands for Leukocyte Ig-like Receptors. J Immunol 2016, 196(3):947-955.
31. Bouloc A, Bagot M, Delaire S, Bensussan A, Boumsell L: Triggering CD101 molecule on human cutaneous dendritic cells inhibits T cell proliferation via IL-10 production. Eur J Immunol 2000, 30(11):3132-3139.
32. Runyan CA, Downey-Slinker ED, Ridpath JF, Hairgrove TB, Sawyer JE, Herring AD: Feed Intake and Weight Changes in Bos indicus-Bos taurus Crossbred Steers Following Bovine Viral Diarrhea Virus Type 1b Challenge Under Production Conditions. Pathogens 2017, 6(4).
33. Afgan E, Baker D, van den Beek M, Blankenberg D, Bouvier D, Cech M, Chilton J, Clements D, Coraor N, Eberhard C et al: The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2016 update. Nucleic Acids Res 2016, 44(W1):W3-W10.
34. Kim D, Langmead B, Salzberg SL: HISAT: a fast spliced aligner with low memory requirements. Nat Methods 2015, 12(4):357-360.
35. Liao Y, Smyth GK, Shi W: featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30(7):923-930.
36. Soneson C, Love MI, Robinson MD: Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res 2015, 4:1521.
37. Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, Haw R, Jassal B, Korninger F, May B et al: The Reactome Pathway Knowledgebase. Nucleic Acids Res 2018, 46(D1):D649-D655.
38. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, Vilo J: g:Profiler-a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res 2016, 44(W1):W83-89.
39. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP et al: STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 2015, 43(Database issue):D447-452.
40. Benjamini Y, Hochberg Y: Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B (Methodological) 1995, 57(1):289-300.

Tables

 Table 1. Differentially expressed genes and their functions identified between vaccine treatment groups during sub-clinical infection with bovine viral diarrhea virus

Table 2. Gene ontology and pathway analysis in sub-clinical animals following BVDV infection

Additional File

Excel workbook showing list of gene expression values (FDR ≤0.15) for each of the comparisons of  KV vs. NON, MLV vs. NON, and MLV vs. KV.