Animals were housed and samples collected in accordance with the Animal Welfare Act Amendments (7 U. S. Code e § 2,131 to § 2,156). All protocols and procedures were approved by the Institutional Animal Care and Use Committee of the National Animal Disease Center, Ames, IA (ARS-22-1035) before the onset of the study. Animal caretakers, veterinarians, and scientific staff performed animal monitoring and record keeping. Animals were humanely euthanized by intravenous sodium pentobarbital following the per-label dose and the discretion of the clinical veterinarian.

The animals were purchased from a commercial facility by USDA. The present study uses tissue and blood samples collected from the same animals in a previous study that has been published (13).

Nine bison (~1–5 years old) were randomly assigned to one of two treatment groups: unvaccinated control (n = 4) or vaccinated (n = 5). An injectable subunit vaccine was formulated by mixing 100 μg of recombinant EFTu and 100 μg of Hsp70 with an oil-in-water emulsified adjuvant (Emulsigen®-D) as previously described (13, 14). Unvaccinated control animals were given a mock vaccine containing only emulsified adjuvant. Two doses of either vaccine or mock vaccine (2 mL per dose) were given at day 0, and at 21 days post-vaccination (DPV). Animals were inoculated with bovine herpesvirus-1 (BHV-1) at 36 DPV and with M. bovis at 40 DPV using an intranasal mucosal atomization device. Inoculums were prepared as previously described (5). Animals were euthanized between 51–65 DPV. Euthanasia and event summary tables can be found in Supplementary Table S1. Peripheral blood was collected by venipuncture and frozen in PAXGene tubes at 5 time points: Day 0 (1st vaccination), 21 DPV (2nd vaccination), 36 DPV (BHV-1 inoculation), 40 DPV (M. bovis inoculation), and 47 DPV (1 week post M. bovis inoculation). Blood samples were taken at the same time for all animals and collected directly before vaccination or inoculation events. Liver, palatine tonsil (PT), retropharyngeal lymph node (RPLN), tracheobronchial lymph node (TBLN), and spleen were collected at necropsy and immediately stored in RNAlater (Thermo Fisher Scientific, Waltham, MA, United States).

Total RNA was extracted and purified using the MagMAX™ mirVana™ Total RNA Isolation Kit (Thermo Fisher Scientific) according to manufacturer’s instructions. Quality and concentration of total RNA samples was evaluated using the Series II RNA 6000 Nano LabChip® kit (Agilent Technologies, Santa Clara, CA, United States) and samples with an RNA integrity number of ≥7.5 were used for subsequent library preparation.

The NEBNext Poly(A) mRNA Magnetic Isolation Module was used for mRNA selection from extracted total RNA and libraries were prepared using the NEBNext Ultra II RNA library prep kit (New England Biolabs, Ipswitch, MA, United States) according to manufacturer’s instructions. Library index incorporation was done using the NEBNext Multiplex Oligos for Illumina kit (New England Biolabs). Libraries were pooled in equal concentration and stored at −20°C until sequencing using the NovaSeq6000 system on S2 flow cells with 150-cycle paired end reads (2 × 150 bp) at the Iowa State University DNA sequencing facility.

Quality of raw reads was assessed with FastQC (version 0.12.1). Adapter sequences were removed and low-quality bases were trimmed (quality-cutoff of 30 and minimum-length ≥60) using Cutadapt (version 4.0). The Bison genome fasta file and annotation files (Bison_UMD1.0) were downloaded from Ensembl. For read alignment and quantification of transcript abundance, the STAR-RSEM pipeline was used (STAR version 2.7.10b; RSEM version 1.3.3). Mapping indices were created and trimmed reads were mapped with STAR. Transcript expression was quantified with rsem-calculate-expression and rsem-generate-data-matrix was used to generate the transcript count matrix.

The software DESeq2 (version 1.44.0) was used for differential expression analysis. The DESeqDataSetFromMatrix and estimateSizeFactors functions were used to create a DESeq2 object and normalize counts with the median ratio method, respectively. Transcripts with at least 5 normalized counts in at least 2 samples were kept for further analysis. The Wald test was used to generate baseMean, log2FoldChange, log fold change standard error, and adjusted p-values. The Benjamini-Hochberg procedure was used by DESeq2 for p-value adjustment in order to control the false discovery rate. Transcripts were deemed significant with an adjusted p-value < 0.05.

Correlation heatmaps for all samples were generated with pheatmap (version 1.0.12). Volcano plots and heatmaps were generated using the R packages EnhancedVolcano (version 1.22.0) and ggplot2 (version 3.5.1). Venn diagrams were generated using ggVennDiagram (version 1.5.2). Gene ontology and KEGG enrichment analysis of DETs was performed with DAVID.

The WGCNA R package (version 1.73) was used to identify modules of highly correlated transcripts connected to vaccination status (15). The top (75% quantile) transcripts for expression variance were identified using the quantile function and used for module identification. The weighted network adjacency matrix was calculated from the expression data based on Pearson correlation coefficients using the adjacency function and a soft-thresholding power of 5 was used to achieve a scale free topology index of 0.85. The adjacency matrix was converted into a topological overlap matrix with the TOMsimilarity function and hierarchical clustering grouped transcripts into modules with a minimum size of 10 transcripts. Eigen values were correlated to vaccination status using cor and corPvalueStudent functions of WGCNA. Network files were exported and visualized using Cytoscape (version 3.10.3).

To quantify transcript expression, raw reads were processed and aligned to the bison reference genome (Bison_UMD1.0). An average of ~36 million raw reads were obtained per sample. Adapter and quality trimming yielded an average of ~27 million trimmed reads per sample. In tissues samples, alignment rates ranged from 84.6% (PT) to 88.82% (RPLN). In blood, alignment rates ranged from 82.7% (2nd vaccination) to 87.5% (1st vaccination).

Liver, spleen, and blood samples separated into distinct clusters with high correlations across samples (R > 0.9), whereas PT, RPLN, and TBLN often clustered together indicating greater similarity in gene expression (Figure 1).

A summary of differentially expressed transcripts between control and vaccinated animals in all tissues and time points is shown in Table 1. The greatest number of differentially expressed transcripts (DETs) within tissues was found in spleen, where 57 DETs were identified. Galactin-3 (LGALS3) was downregulated and Fos Proto-Oncogene, AP-1 Transcription Factor (FOS) was upregulated in spleen of vaccinated animals (Figure 2). The 2nd highest number of DETs in tissue was found in liver, where major facilitator super family domain containing 2a (MFSD2A) was upregulated in vaccinated animals compared to control (16). The number of DETs in other tissues ranged from 3 (RPLN) to 12 (TBLN). DETs included SEC16 homolog A (SEC16A), which was upregulated in TBLN of vaccinated animals (17). Upregulation of general control non-derepressible protein 1 (GCN1) was observed in RPLN of vaccinated animals. A regulator of interferon signaling, zinc finger E-box binding homeobox 1 (ZEB1), was downregulated and an immunoregulatory molecule, Disabled-2 (DAB2), was upregulated in PT of vaccinated animals (18–20).

The number of DETs between control vs. vaccinated animals in blood ranged from 11 total DETs (Day 0) to 123 total DETs (36 DPV) (Table 1). The transcription factor paired box 5 (PAX5) facilitates differentiation of B cells and was upregulated in vaccinated animals at 21 DPV as shown in Figure 3 (21). Transcripts related to heparin binding, such as apolipoprotein E (APOE) and fibronectin 1 (FN1) were downregulated in vaccinated animals at 36 DPV (22, 23). At 40 DPV, v-set immunoregulatory receptor (VSIR) was upregulated in vaccinated animals (24). Anoctamin 10 (ANO10) was downregulated in vaccinated animals compared to control (25). Sorcin (SRI) and Fc fragment of IgE receptor Ig (FCER1G) were downregulated at 1 week post M. bovis infection (47 DPV) in vaccinated animals (26, 27). DESeq2 results for all samples can be found in Supplementary Tables S2, S3.

Out of 190 DETs, 178 were unique to day of vaccination in blood. However, LOC105004569 (ortholog of NK-lysin) and LOC105001075 (ortholog of T cell-interacting activating receptor on myeloid cells 1/TARM1) were shared across more than one time point (Figure 4, left). The antimicrobial peptide, NK-lysin, is produced by T lymphocytes and natural killer cells, and was upregulated in vaccinated animals at 21 and 36 DPV (28). The leukocyte receptor, TARM1, is critical for macrophage activation and was upregulated in vaccinated animals at day 0 and 36 DPV (29).

In tissue comparisons, DETs in liver overlapped with DETs in TBLN, RPLN, and PT (Figure 4, right). Granulin precursor (GRN) and ribosomal protein L14 (RPL14) were upregulated in liver and TBLN of vaccinated animals. The gene GCN1 was downregulated in liver and upregulated in RPLN of vaccinated animals. Nuclear receptor subfamily 1 group D member 2 (NR1D2) was downregulated in liver and PT of vaccinated animals (30). A negative regulator of Toll-like receptor signaling, DAB2, was upregulated in PT and downregulated in TBLN of vaccinated animals (31).

As shown in Table 2, significant enrichment of DETs was observed in spleen and liver. No significant enrichment was observed in PT, TBLN, or RPLN. DETs in spleen, such as T-complex 1 (TCP1) and prostaglandin E synthase 3 (PTGES3), were downregulated in vaccinated animals and enriched in biological processes of protein stabilization and positive regulation of telomere maintenance. There was a tendency for enriched innate immune response (p-value = 0.06), where LGALS3 and LOC105003308 (Ortholog of interferon-induced guanylate-binding protein/GBP-1) were downregulated in vaccinate spleen and ADAM metallopeptidase domain 15 (ADAM15) was upregulated. In liver, protein tyrosine phosphatase non-receptor type 6 (PTPN6) and beta-transducin repeat containing E3 ubiquitin protein ligase (BTRC) were upregulated in vaccinated animals and enriched in the negative regulation of T cell receptor signaling pathway.

The greatest number of significantly enriched gene ontology terms and pathways was found at 36 DPV (Figure 5). The most significantly enriched pathways were cell adhesion molecules, Th17 cell differentiation, antigen processing and presentation, and HIV-1 infection pathways. Infiltrated lymphocyte markers (GATA binding protein 3/GATA3, Interleukin 27 receptor subunit alpha/IL27RA, CD3E, CD4, and T-box transcription factor 21/TBX21) were enriched in Th17 cell differentiation and upregulated in vaccinated animals at 36 DPV (Figure 6A). The ortholog of bovine leukocyte antigens class I major histocompatibility complex (MHC), alpha chain BL3-7 (LOC104983816) and a binding partner of MHC Class I molecules, TAP binding protein (TAPBP), were upregulated in vaccinated animals and enriched in antigen processing and presentation (Figure 6B). Phospholipase C, gamma-2 (PLCG2) was downregulated and tumor necrosis factor receptor superfamily member 1B (TNFRSF1B) was upregulated in vaccinated animals and enriched in HIV-1 infection (Figure 6C). Upregulation of cell adhesion DETs was found, such as CD2, CD6, LOC105004586 (Ortholog of CD8b), and LOC104986653 (Ortholog of Sialophorin/SPN) at 36 DPV (Figure 6D).

At 40 DPV (or 4 days post-BHV-1 challenge), VPS50 subunit of EARP (VPS50) and SEC31 Homolog A (SEC31A) were upregulated in vaccinated animals and enriched in protein transport. In addition, two upregulated (Fibulin 5/FBLN5 and SPARC Like 1/SPARCL1) and two downregulated (RAS Guanyl Releasing Protein 4/RASGRP4 and S100 Calcium Binding Protein A4/S100A4) DETs in vaccinated animals were enriched in calcium ion binding at 40 DPV.

Weighted co-expression analysis identified 5 modules of co-expressed transcripts with correlations to treatment group at 40 and 47 DPV. Of these modules, Module V displayed a significant (p-value < 0.05) and intermediate, negative correlation (R = −0.49) with vaccinated animals (Figure 7A). Extracellular matrix protein 1 (ECM1), vimentin (VIM), and calpain small subunit 1 (CAPNS1) were among the transcripts with the highest significance (R > 0.7) for vaccination status and connectivity to other transcripts in the module (Figure 7B). The transcripts with the highest overall module membership in Module V included G protein subunit beta (GNB2), Src Family Tyrosine Kinases (HCK and FGR), solute carrier family 15 member 3 (SLC15A3), integrin subunit beta 7 (ITGB7), and flotillin 1 (FLOT1) as shown in Figure 7C. The most significant KEGG pathways enriched with the co-expressed transcripts in Module V included osteoclast differentiation, neutrophil extracellular trap formation, and chemokine signaling pathway (Figure 7D). The most significant enrichment for biological processes and molecular functions included protein localization to plasma membrane, cytokine-mediated signaling pathway, innate immune response, guanyl-nucleotide exchange factor activity, and cytokine receptor activity (Figure 7D).