Animal resources for the present study were obtained from a 26-week vaccination and challenge experiment of age- and sex-matched cattle infected with M. bovis (40–44). Ten male age-matched Holstein–Friesian calves (4–6 months old) were sourced from farms known to be free of BTB disease. The animals used for the experimental work described here were the naive control group (non-vaccinated) for a vaccine efficacy study (40). Figure 1 shows the experimental schedule used by Dean et al. (40) and details the sampling time points for the 10 non-vaccinated control cattle used for the research work described here.

The challenge strain, M. bovis AF2122/97 (2, 45), was delivered endobronchially at 2 × 10³ colony-forming units per animal using the following procedure described by Whelan et al. (46). Prior to endobronchial inoculation animals were sedated with Rompun® (Bayer Animal Health, Newbury, UK) according to the manufacturer's instructions. Following this, an LSVP 22 VGS89x14 endoscope (Veterinary Endoscopy Services, Welshpool, UK) lubricated with Vet-Lubigel (Millpledge Veterinary, Clarborough, UK) was inserted through a nostril into the trachea and placed above the bronchial opening to the cardiac lobe and the main bifurcation between left and right lobes. A cannula of 1.8-mm internal diameter (Veterinary Endoscopy Services) was inserted through the endoscope and used to deliver the M. bovis AF2122/97 inoculum in 2 mL of phosphate-buffered saline (PBS). Following this, an additional 2 mL of PBS was then used to wash the cannula, and the cannula and endoscope were withdrawn. For each individual animal, a new sterile cannula was used, the internal channel of the endoscope, through which the cannula was inserted, was rinsed with 20 mL of PBS, and the outside of the endoscope was cleaned thoroughly with sterilising tissue wipes (Medichem International, Sevenoaks, UK).

Individual responses to infection across the time course and disease pathology for the animals used in this study have been described in detail previously and include whole-blood IFN-γ assay, evaluation of peripheral blood mononuclear cell (PBMC) cytokine responses by intracellular cytokine staining, gross (visible) pathology and histopathology, and evaluation of bacterial load in lymph nodes (40, 41).

Approximately 3 mL of ex vivo peripheral blood was sampled from all 10 naive control animals at −1 week pre-infection and then at +1, +2, +6, +10, and +12 weeks post-infection (Figure 1). All blood samples were obtained during the morning (between 7:00 and 10:00 a.m.) of each collection day and directly collected into Tempus™ blood RNA tubes (Applied Biosystems®/Thermo Fisher Scientific, Warrington, UK). Immediately after blood collection at each time point, Tempus™ tube samples for each animal were vortexed for ~10 s to ensure complete red blood cell lysis. Tempus™ tube blood lysate samples for animals at each of the nine time points were then stored at −80°C until they were used for total RNA extraction and purification.

The Tempus™ Spin RNA Isolation Kit (Applied Biosystems®/Thermo Fisher Scientific) was used for total RNA extraction and purification using the following protocol provided by the manufacturer. Tempus™ tube blood lysate samples were thawed at room temperature prior to RNA extraction and purification. Once thawed, for each sample, ~3 mL of blood lysate was transferred to a 50-mL plastic centrifuge tube, and PBS was added to a final volume of 12 mL. Each sample was then mixed by vortexing for 30 s and then centrifuged at 3,000 × g for 30 min at 4°C. The supernatant was then removed, and the remaining RNA-containing pellet was resuspended with a brief vortex in 400 μL of the proprietary RNA Purification Resuspension Solution. Following this, the resuspended RNA sample was pipetted into the RNA purification filter inserted into a 1.5-mL microcentrifuge tube for waste collection. The RNA purification filter/microcentrifuge tube was then centrifuged at 16,000 × g for 30 s, and the liquid waste and microcentrifuge tube discarded. The RNA purification filter was then placed in a clean microcentrifuge tube, 500 μL of proprietary RNA Purification Wash Solution 1 was added, followed by another centrifugation step at 16,000 × g for 30 s and disposal of the liquid waste and microcentrifuge tube. This step was then repeated using 500 μL of proprietary RNA Purification Wash Solution 2 with a centrifugation step at 16,000 × g for 30 s. A final wash step was then performed with 500 μL of RNA Purification Wash Solution 2 and centrifugation at 16,000 × g for 30 s followed by disposal of the liquid waste and microcentrifuge tube. The RNA purification filter was then placed in a clean microcentrifuge tube and centrifuged at 16,000 × g for 30 s to dry the membrane. The RNA purification filter was then inserted into a clean RNase-free collection microcentrifuge tube and 100 μL of Nucleic Acid Purification Elution Solution was added and incubated for 2 min followed by centrifugation at 16,000 × g for 30 s; the RNA eluate was then pipetted back onto the filter membrane, and the centrifugation step was repeated. Approximately 90 μL of the final RNA eluate was then pipetted (avoiding particulate material) into a new labelled RNase-free collection microcentrifuge for long-term storage at −80°C.

RNA quantity and quality checking were performed using a NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer using an RNA 6000 Nano LabChip kit (Agilent Technologies, Cork, Ireland). The majority of samples displayed a 260/280 ratio >1.8 and RNA integrity numbers >8.0 (Supplementary Table 1 in Supplementary Material 1). RNA quality and quantity checking revealed that three samples did not have measurable quantities of RNA, and these were excluded from downstream RNA-seq library preparation (15, ID 6520, +2 weeks; 21, ID 6522, +2 weeks; and 27, ID 6526, + 2 weeks).

For RNA-seq library preparation, 1 μg of total RNA from each sample was used to prepare individually barcoded strand-specific RNA-seq libraries. Two rounds of poly(A)⁺ RNA purification were performed for all RNA samples using the Dynabeads® mRNA DIRECT™ Micro Kit (Ambion®/Thermo Fisher Scientific, Austin, TX, USA) according to the manufacturer's instructions. The purified poly(A)⁺ RNA was then used to generate strand-specific RNA-seq libraries using the ScriptSeq™ v2 RNA-Seq Library Preparation Kit, the ScriptSeq™ Index PCR Primers (sets 1–4), and the FailSafe™ PCR enzyme system (all sourced from Epicentre®/Illumina® Inc., Madison, WI, USA) according to the manufacturer's instructions. RNA-seq libraries were purified using the Agencourt® AMPure® XP system (Beckman Coulter Genomics, Danvers, MA, USA) according to the manufacturer's instructions for double size selection (0.75 × followed by 1.0 × ratio). RNA-seq libraries were quantified using a Qubit® fluorometer and Qubit® dsDNA HS Assay Kit (Invitrogen™/Thermo Fisher Scientific, Carlsbad, CA, USA), whereas library quality cheques were performed using an Agilent 2100 Bioanalyzer and High Sensitivity DNA Kit (Agilent Technologies Ltd.). Individually barcoded RNA-seq libraries were pooled in equimolar quantities, and the quantity and quality of the final pooled libraries (three pools in total) were assessed as described previously. RNA-seq library sample barcode index sequences are detailed in Supplementary Table 1 (Supplementary Material 1).

Prior to high-throughput sequencing, the content of several RNA-seq libraries was validated using conventional Sanger dideoxy sequencing. Library inserts from 16 libraries were cloned using the Zero Blunt® TOPO® PCR Cloning Kit according to the manufacturer's instructions (Invitrogen™/Thermo Fisher Scientific). Sanger sequencing of 36 plasmid inserts from these selected libraries confirmed that the RNA-seq libraries contained inserts derived from bovine mRNA. Plasmid sequencing was outsourced (Source Bioscience Ltd., Dublin, Ireland), and sequences generated were validated using BLAST searching of the DNA sequence database (47). Cluster generation and high-throughput sequencing of the pooled RNA-seq libraries were performed using an Illumina® HiSeq™ 2000 Sequencing System at the MSU Research Technology Support Facility (RTSF) Genomics Core (https://rtsf.natsci.msu.edu/genomics; Michigan State University, MI, USA). Each of the three pooled libraries was sequenced independently on five lanes split across multiple Illumina® flow cells. The pooled libraries were sequenced as paired-end 2 × 100 nucleotide reads using Illumina® version 5.0 sequencing kits. Additionally, after exploratory data analysis (Supplementary Figures 1, 2), it was decided to remove animal ID 6522 completely from the analysis and proceed with 52 RNA-seq sample data sets (Supplementary Table 1 in Supplementary Material 1). All RNA-seq data generated for this study have been deposited in the European Nucleotide database with experiment series accession numbers (PRJEB27764 and PRJEB44470).

Except where indicated, bioinformatics procedures and analyses were performed on a 32-core Compute Server running Linux Ubuntu (version 12.04.2) hosted at the UCD Research IT Data Centre (stampede.ucd.ie) and administered by the UCD Animal Genomics Group. All of the bioinformatics workflow/pipeline components including Linux Bash, Perl, and R scripts used were deposited in a GitHub repository (https://github.com/kmcloughlin1/RNA-sequencing) and were modified from published methods described by our group (48). Supplementary Figure 3 shows a schematic of the complete RNA-seq bioinformatics workflow and the downstream tools used for time-series analysis and various systems biology methods.

Deconvolution (filtering and segregation of sequence reads based on the unique RNA-seq library barcode index sequences; Supplementary Table 1 in Supplementary Material 1) was performed by the MSU RTSF Genomics Core using a pipeline that simultaneously demultiplexed and converted pooled sequence reads to discrete FASTQ files for each RNA-seq sample with no barcode index mismatches permitted. The RNA-seq FASTQ sequence read data were then downloaded from the MSU RTSF Genomics Core FTP server, and a custom Perl script was used to filter out paired-end reads containing adapter sequence contamination (with up to three mismatches allowed) and to remove poor quality paired-end reads (i.e., one or both reads containing 25% of bases with a Phred quality score <20). The quality of individual RNA-seq sample library files was then reassessed postfiltering using the FastQC software package (version 0.10.1) (49).

Paired-end reads, from each filtered individual library, were aligned to the Bos taurus reference genome (UMD3.1.73) (50) using the STAR aligner software package (version 2.3.0) (51). For each library, raw counts for each gene based on the sense strand data were obtained using the featureCounts software from the Subread package (version 1.3.5-p4) (52). The featureCounts parameters were set to unambiguously assign uniquely aligned paired-end reads in a stranded manner to the exons of genes within the UMD3.1.73 B. taurus reference genome annotation. The gene count outputs were then used to perform differential gene expression analysis using the edgeR Bioconductor package (version 3.2.4) (53) within an R-based pipeline that was customised to perform the following functions:

The Ingenuity® Pathway Analysis (IPA) software package (57) with the Ingenuity® Knowledge Base (Qiagen Corp., Redwood City, CA, USA; release date July 2014; www.ingenuity.com) was used to identify overrepresented (enriched) canonical pathways and construction of biological interaction networks for sets of DE genes at each post-infection time point (+1, +2, +6, +10, and +12 weeks) compared to the pre-infection time point (−1 week). For identification of overrepresented canonical pathways, a multiple testing correction (Benjamini–Hochberg method) was applied with an FDR threshold ≤ 0.05. The IPA Biomarker Filter tool was also used to identify and prioritise molecular biomarker candidates such that only experimentally observed and high-confidence predicted biological relationships were included. All the IPA data sources were used for three mammalian species in the IPA Knowledge Base (Homo sapiens, Mus musculus, Rattus norvegicus). Biological interaction networks were ranked according to Network Score values generated algorithmically by IPA and based on the hypergeometric distribution and calculated with the right-tailed Fisher exact test (58).

Time-series analysis of gene expression data from the animal infection time-course experiment was performed using the Short Time-series Expression Miner (STEM) software package (59). The computational procedure for selecting model profiles that are representative and distinct is described by Ernst et al. (60). The software package implements a method for clustering short time-series expression data that can differentiate between real and random patterns of temporal gene expression changes and assigns each gene to the model profile that most closely matches the temporal gene expression profile for that gene as determined by the correlation coefficient. A permutation test is then used to determine which model profiles have a statistically significant number of genes assigned compared to random expectations from the mean number assigned to each profile based on the permuted data (59). STEM also incorporates GO enrichment functionality for biological interpretation of time-series gene expression data.

Deconvolution and filtering of sequence reads to remove adaptor-dimer contamination yielded a mean of 20.6 ± 2.0 million reads per individual barcoded RNA-seq sample library (n = 52 libraries and ± SD); this corresponded to a mean of 82% reads that passed this filtering step. These filtered reads were then aligned to the B. taurus UMD3.1.73 genome build. This yielded a mean of 15.4 ± 1.7 million filtered reads (91.5%) that uniquely mapped to this bovine genome build with a mean mapped length of 195.6 ± 0.6 bp; a mean of 0.77 ± 0.19 million reads (4.6%) that mapped to multiple genomic locations and 0.67 ± 0.17 million reads (3.9%) that did not map to any genomic location. Further analysis demonstrated that a mean of 63.1% of the filtered uniquely mapping reads (9.7 ± 1.1 million reads) were assigned to Ensembl gene IDs for the UMD3.1.73 genome build and 36.9% (5.7 ± 0.86 million reads) were ambiguous or could not be assigned to an annotated genomic region. Supplementary Table 2 (Supplementary Material 1) and Supplementary Figure 4 show the RNA-seq summary statistics. Filtering of the RNA-seq data using 52 samples (60 minus the three technical dropouts and the animal ID 6522 samples at −1, +1, +6, +10, and +12 weeks) produced 12,406 genes suitable for downstream differential expression analysis.

Multidimensional scaling plots (Supplementary Figure 5) demonstrated that it was not possible to differentiate the infected animals from non-infected animals during the early stages of the animal infection time course (+1 week post-infection vs. −1 week pre-infection, and +2 weeks vs. −1 week pre-infection). Conversely, discrimination of infected and non-infected animals was partially observed at +6 weeks post-infection and was clearly evident at +10 weeks post-infection, but this pattern of discrimination effectively disappeared by +12 weeks post-infection. Previous work by our group has shown that microarray and RNA-seq gene expression data sets from peripheral blood leukocytes (PBLs) can be used to unambiguously discriminate M. bovis–infected from non-infected control cattle (28, 31). However, it is important to note that the M. bovis–infected animals used for these earlier studies were heavily infected animals, which were maintained for ongoing disease surveillance and potency testing of diagnostic reagents.

The concept and the implications regarding biomarker identification and biosignature development for infectious disease have been explored thoroughly by Chaussabel et al. (36, 61, 62). In particular, these researchers emphasise that leukocytes present in peripheral blood convey valuable information about the status of the immune system that can be translated to biomarkers and onward to sensitive and specific biosignatures of infection. In addition, peripheral blood is easily accessible and can be stabilised and processed for high-throughput transcriptomics analyses using RNA-seq and other massively parallel gene expression technologies such as microarrays. It is also notable that development of transcriptional biomarkers in cattle is relatively straightforward because very low levels of globin gene transcripts (HBA and HBB) are observed in bovine peripheral blood compared to other mammalian species (63, 64).

Statistical analysis of the RNA-seq gene expression data with a B-H FDR adjusted P ≤ 0.05 demonstrated that differential gene expression was evident between each of the five post-infection time points (+1, +2, +6, +10, and +12 weeks) and the −1 week pre-infection time point (Figure 2A). Relatively small numbers of DE genes were detected at +1 week (37 exhibited increased and 20 exhibited decreased expression) and +2 weeks (83 increased and 10 decreased); however, the numbers of DE genes were substantially greater at +6 weeks (415 increased and 272 decreased), +10 weeks (1,278 increased and 1,305 decreased), and +12 weeks (222 increased and 116 decreased). Supplementary Tables 3–7 (Supplementary Material 2) provide detailed information on the differential expression analysis results at each of the five post-infection time points relative to the −1 week pre-infection time point for the complete set of 12,406 filtered genes. Figure 2B also shows a Venn diagram for the significant DE genes at each of the five post-infection time points (relative to the −1 week pre-infection control time point). In addition, Table 1 provides detailed information for 19 genes shown in Figure 3 that were significantly DE across all the five post-infection time points compared to the −1 week pre-infection control time point.

As shown in Figure 3, there is a striking concordance between the patterns of expression for these 19 genes across the infection time course and earlier PBL microarray and RNA-seq studies published by our group (28, 31). These results provide good support for the hypothesis that a biosignature of M. bovis infection can be generated using transcriptomics data from cattle with early- and later-stage BTB. It also provides evidence that putative transcriptional biomarkers, identified using an experimental challenge with a relatively high M. bovis infectious dose, can be translated as diagnostic tools for use in naturally infected animals. Diagnostic biosignature development focusing on smaller panels of transcriptomics-based biomarkers has been used with notable success for human TB. In this case, research work has focused on specificity and differentiating active TB from latent TB and also TB disease from non-infected controls and diseases with similar pathology but distinct aetiology such as sarcoidosis, pneumonia, and lung cancer (35, 39, 65–71).

Previous work on M. tuberculosis and human immunodeficiency virus infection in humans has shown that CXCR4 is upregulated in blood monocytes and bronchoalveolar lavage cells from human patients with pulmonary TB (72–74). In addition, a loss-of-function mutation in the murine Thbd gene that impairs activated protein C production results in uncontrolled lung inflammation in mice infected with M. tuberculosis, highlighting the importance of the THBD gene in mammalian TB disease (75). Galietti et al. have also shown that M. tuberculosis– and M. bovis–infected, but not M. avium–infected, human monocytes showed increased expression of the CDKN1A protein encoded by CDKN1A (76). A range of studies have shown that levels of the protein product of the PLAUR gene are elevated in serum from human patients infected with M. tuberculosis (77–81). Increased expression of the OSM gene and induction of matrix metalloproteinases, which contribute to tissue damage characteristic of TB, have been demonstrated in M. tuberculosis–infected human monocytes (82). High expression of OSM was also observed in the blood transcriptome of patients presenting with high mycobacterial load sputum (83). In addition, the OSM gene is located within a candidate QTL region for TB susceptibility identified using admixture mapping in humans (84).

Ten of the 19 genes that showed consistently increased expression across all post-infection time points were observed to overlap with results from RNA-seq of an in vitro infection time-course experiment using bovine AMs stimulated with M. bovis (48). FOSB and NR4A1 were upregulated in AMs at 2 hpi; CXCL8, NR4A1, PLAUR, and RGS16 were upregulated at 6 hpi; EVI2A, CXCL8, FOSB, HBEGF, OSM, PLAUR, RGS16, and THBD were upregulated at 24 hpi; the eight genes observed at 24 hpi plus the CDKN1A gene were also upregulated at 48 hpi.

One of the most notable putative transcriptional biomarkers represented in the panel of 19 genes and in the independent PBL studies is CXCL8 (previously known as IL8). CXCL8 is a chemokine encoded by the CXCL8 gene, which is a strong neutrophil chemoattractant and also chemotactic for monocytes and T cells (85, 86); it has been observed to exhibit increased expression for different mycobacterial infections in a range of mammalian systems (87–95). CXCL8 enhances killing of mycobacteria by neutrophils and macrophages (96, 97), and these immune cells also secrete CXCL8 when stimulated by M. tuberculosis (98). In this regard, Godaly and Young showed that M. bovis bacillus Calmette–Guérin (BCG) induces CXCL8 secretion by human neutrophils via MyD88-dependent signalling through TLR2 and TLR4 (99). Also, stimulation of human lung fibroblasts in vitro using conditioned medium from M. tuberculosis–infected monocytes caused prolonged expression of CXCL8 mRNA and >10-fold increase in CXCL8 secretion (88).

With regard to the CXCL8 mRNA transcript as a biomarker of infection, Alessandri et al. were able to detect significantly elevated levels of the CXCL8 cytokine in plasma from patients with pulmonary TB (87). More recently, based on reversion of IGRA test results in a Chinese cohort, it has been proposed that decreased serum levels of CXCL8 are associated with clearance of M. tuberculosis infection (100). In addition, using microarray and reverse transcriptase–qualitative polymerase chain reaction technologies, Widdison et al. have shown that M. tuberculosis– and M. bovis–infected bovine AMs express high levels of CXCL8 transcripts compared to non-infected control cells (90). In support of this, using RNA-seq, we have shown that CXCL8 increases in expression in bovine AMs infected with either M. tuberculosis or M. bovis across a 48-h time course (101). CXCL8 has also been shown to be significantly increased in expression after in vitro PPD-b stimulation of PBMCs from cattle infected with M. bovis (102) and bovine monocyte-derived macrophages (103). CXCL8 also exhibited increased expression in PBL (28) but decreased expression in non-stimulated PBMCs from M. bovis–infected cattle (104). Also, Almeida de Souza et al. have shown that antimycobacterial treatment reduces high plasma levels of CXCL8 and other CXC chemokines detected in plasma from human patients with active TB (91), and Huang et al. have also demonstrated that AMs and PBMCs from TB patients express CXCL8 at significantly higher levels than healthy controls (94). Interestingly, the potential specificity of increased CXCL8 gene expression as a biomarker for M. bovis infection in cattle is illustrated by recent results obtained by Alonso-Hearn et al. (105). Using similar RNA-seq methodology, they observed significantly decreased expression of CXCL8 in peripheral blood from cattle infected with M. avium subsp. paratuberculosis, the causative agent of Johne disease. Finally, it is important to note that several primary studies and meta-analyses have provided evidence that a single-nucleotide polymorphism (rs4073) at the human CXCL8 gene locus is associated with resistance/susceptibility to M. tuberculosis infection (106–109).

Of the 12,406 genes (50.40% of total B. taurus reference genes) that were suitable for differential expression analysis, 10,703 genes (86.27%) were mapped to molecules in the IPA Knowledge Base. IPA was used to identify overrepresented canonical pathways and construct biological interaction networks for sets of DE genes at each post-infection time point (+1, +2, +6, +10, and +12 weeks) compared to the pre-infection control time point (−1 week). Only DE genes that were significant after a multiple testing correction was applied (Benjamini–Hochberg method, FDR threshold ≤ 0.05) were used. The gene expression data for the panel of 19 genes (Table 1) were also analysed using IPA and to identify enriched canonical pathways and biological interaction networks. However, these analyses did not reveal any notable functionally relevant pathways or networks (results not shown).

The STEM tool was designed specifically for analyses of short time-series data sets (3–8 time points) (59, 60) similar to the RNA-gene expression data set obtained from the M. bovis animal infection time-course experiment described here. Time-series analysis can be a powerful technique for uncovering networks of coregulated genes in longitudinal time-course experiments (136–138), particularly for gene expression data associated with host immunobiological responses to infection (139–142).

For the present study, STEM time-series analyses of differential gene expression across the M. bovis infection time course demonstrated that large groups of genes exhibited comparable patterns of gene expression across the five post-infection time points (+1, +2, +6, +10, and +12 weeks post-infection) relative to the −1 week pre-infection time point. Two different STEM analyses were performed based on (a) expression data for all detectable expressed genes across the infection time course, which corresponded to 4,103 genes (STEM analysis 1) and (b) the union set of DE genes from all post-infection time points vs. −1 week pre-infection, which corresponded to 2,935 genes (STEM analysis 2). Supplementary Figure 8 shows the top 50 time-series model profiles obtained for STEM analysis 1 and STEM analysis 2.

The top-ranked time-series profiles (by P-value) obtained using the two different STEM analyses were very similar. For example, as shown in Supplementary Table 12 (Supplementary Material 4), the first-ranked STEM model profile for STEM analysis 1 (profile 40; Figure 7) was enriched for the signal transduction, single organism signalling, cell communication, regulation of multicellular organismal process, and cellular response to stimulus GO terms. This profile was also similar to the second-ranked model profile for STEM analysis 2 (profile 40; Supplementary Figure 9), which was enriched for many of the same GO terms (signal transduction, single organism signalling, cell communication; see Supplementary Table 13, Supplementary Material 4). The third-ranked model profile for STEM analysis 1 (profile 23; Supplementary Figure 9) was highly similar to the first-ranked model profile for STEM analysis 2 (profile 23; Figure 7) with exactly the same overrepresented GO terms (mitochondrial inner membrane, organelle inner membrane, mitochondrial part, mitochondrial membrane, mitochondrial envelope; Supplementary Tables 14, 15, Supplementary Material 4).

It is interesting to note that STEM profile 40 in each analysis is characterised by a cluster of ~300 genes associated with cell signalling and cellular response to stimuli, which exhibited increasing expression across the four time points with a peak at +10 weeks, followed by a substantial decrease at +12 weeks (Figure 7). Conversely, STEM profile 23 (516 genes for STEM analysis 2) is characterised by genes associated with mitochondrial components, particularly the mitochondrial membrane, which displayed an oscillating pattern of expression with a marked decrease at +1 week 10 post-infection (Figure 7).

These time-dependent patterns of gene expression in peripheral blood may reflect pathogenesis of early BTB disease during the infection time course with concomitant host cellular responses to M. bovis infection, disruption of homeostasis, and changing cellular, tissue, and organismal energy requirements (143–145). In addition, it is important to note that although these longitudinal patterns of gene expression may be due to coregulation of genes in the same cluster, they are likely to also reflect fluctuations in peripheral blood cell type populations comparable to those previously observed for comparisons of M. bovis–infected and control non-infected cattle (28, 104).