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

Background: Bovine viral diarrhea (BVD) is a widespread immunosuppressive disease of cattle. Infection with bovine viral diarrhea virus (BVDV) can also be sub-clinical and persistent in herds, leading to high economic impact to cattle producers. Complicating this is BVD’s ability to predispose cattle to bovine respiratory disease (BRD) complex. Although many producers vaccinate for BRD to help protect herds from this predisposition, there are gaps in the knowledge of differences between the protection afforded by modified live (MLV) or killed (KV) vaccines, and the level of cross-protection they offer to BVDV. These knowledge gaps affect the economic impact of BVDV because, with more information on vaccine efficacy, cattle producers could possibly reduce treatment costs and disease spread. The study objective was to determine the level of BRD vaccine cross-protection by examining the underlying differential gene expression response to a sub-clinical BVDV challenge in Nellore-Angus cattle. Results: Statistical analysis of clinical data showed that MLV had lower BVDVa, BVDV2, and BVDVb antibody titers (p < 0.001) compared to the KV vaccine. The three treatment groups (no vaccine (NON), KV, MLV; n = 5 per group) were also compared by RNA-seq analysis of buffy coats from cattle at 14 dpi. The BRD KV group displayed a higher antibody (AB) response, however transcriptomic analysis showed that the KV group still showed high innate gene expression at 14 dpi. The MLV showed displayed a higher number of adaptive and cell-mediated immune response genes being differentially expressed in reaction to challenge at 14 dpi. Analysis also showed the difference in vaccine cross-protection involved evidence of complement activation as part of the MLV bolstered immune response. Lack of detection of complement activation in the KV group could explain the higher titers, while complement activation in the MLV group could lead to less AB response due to better cell-mediated response. Upregulation of MHC molecules and T-cell related gene expression was also observed within the MLV group

only. Conclusions: Although both vaccines provide some level of cross-protection over non-vaccinated animals, the results point to MLV producing cross-protection that extends the adaptive immune response to BVDV.

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][4][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.

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)

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) [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 crossprotection 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 crossprotective 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]  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 proinflammatory 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. 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 immunoglobulinlike 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.

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

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][38][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 ® (

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.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to