Trained individuals within conventional beef production systems handle cattle at a limited number of time points, predominately at arrival, to reduce animal stress and minimize labor demand. Coupled with the complex nature of BRD and lack of a gold standard diagnostic test, these systems have difficulty determining likelihood of disease acquisition and severity on an intra-population scale. Accordingly, our previous RNA-Seq studies were conducted to better characterize the at-arrival host response in cattle related to subsequent development and treatment of naturally-acquired BRD within the first 28 days of facility placement [30, 31, 32]. RNA-Seq studies employ robust statistical and mathematical models in an effort to identify potential molecular targets and mechanisms that may improve disease understanding, detection, and therapy [43, 45]. Such studies reduce the overall dimensionality of host expression and allow for the determination of a finite number of molecules to assess in future studies. However, the development of disease prediction or prognostic models requires rigorous testing and can be resource intensive, as gene-by-environment interactions and disease temporality are difficult to account for in developing predictive models [46, 47]. Therefore, we focused on corroborating differential expression and retrospectively classifying the predictive capability of mRNA molecules previously identified by RNA-Seq, but through a method feasible for use in large numbers of cattle. To this purpose, we utilized NanoString nCounter mRNA profiling as a less variable and higher throughput method of evaluating specific gene expression, compared to the commonly utilized alternative, RT-qPCR [48, 49, 50, 51, 52]. To our knowledge, this study is the first to substantiate selective at-arrival gene expressional patterns previously associated with BRD acquisition and severity, utilizing blood samples from beef cattle across multiple independent populations.
A limitation of this work is that we aimed to identify at-arrival gene expression profiles that predicted future treatment for BRD based on clinical assessment, which can lack diagnostic sensitivity [13, 14]. While we recognize that the visual identification of clinical BRD is relatively insensitive and potentially confounded by study location, our statistical analyses and modeling were conducted to account for this effect. Furthermore, analysis of weight gain records at time of sale (ADGend) demonstrated an objective loss in production that was associated with increased frequency of treatment, in agreement with previous research (Fig. 3) [53, 54], confirming the significance of the BRD diagnoses in this study.
In total, 30 unique genes were determined to be differentially expressed across all comparisons (Fig. 1; Supplemental Tables S4-S7). Cattle never diagnosed with clinical BRD possessed increased expression of genes broadly associated with leukocyte/granulocyte development and differentiation (GCSAML, IL5RA, KLF17, LOC100297044 (CCL14), and LOC100335828 (CD200R1)), SPMs (ALOX15 and HPGD), anti-inflammatory/apoptotic processes (CPB2 and ITGA4), and antimicrobial peptides (CATHL3, GZMB, and LTF). CD200R1, CCL14, GCSAML, and IL5RA have previously been identified as chemokines for macrophage activity and enhanced lymphocyte migration and survival [55, 56, 57, 58, 59]. SPMs, while not widely explored in ungulates, are categorized into two primary categories based on the metabolism of arachidonic acid (lipoxins) or essential polyunsaturated fatty acids (resolvins, maresins, and protectins) [60, 61]. In other mammalian species, SPMs are involved in quelling prolonged pro-inflammatory events and promoting cellular/molecular homeostasis in the lower airways; these molecules are currently under investigation as therapeutic modalities in humans with sepsis or viral-induced pneumonia/acute respiratory distress syndrome (ARDS) [62, 63, 64, 65]. CPB2 and ITGA4 are key regulators of pro-inflammation, as they are involved in modulating the complement system, degrading induced plasma anaphylatoxins, and reducing nitric oxide accumulation [66, 67, 68, 69, 70, 71]. Cattle are capable of producing a myriad of well conserved innate peptides within the cathelicidin and defensin peptide families, such as CATHL3, GZMB, and LTF, which have oxygen-independent killing capacity against both Gram-positive and Gram-negative pathogens [72, 73]. These amphophilic molecules align specifically to bacterial cell walls via electrostatic interactions and penetrate into the cytoplasmic space, leading to DNA replication disruption and/or bacterial autolysis [74, 75, 76, 77]. Collectively, these findings indicate that cattle that remain clinically healthy possessed enhanced immunological, anti-inflammatory, and microbial killing mechanisms at arrival, compared to cattle that require antimicrobial therapy within the first 28 days of arrival.
Cattle that would go on to develop BRD possessed increase expression of genes at arrival broadly associated with pro-inflammatory, granulocytic processes (LRG1, MCF2L, SPP1), alternative complement (CFB), and type I interferon signaling (HERC6, IFI6, ISG15, and MX1), when compared to healthy cattle. Notably, Treated_1 cattle only possessed increased expression of CFB and LRG1 when compared to Healthy cattle. Therefore, these aforementioned genes and their associations are predominantly detected within Treated_2 + cattle. LRG1 encodes for a leucine-rich glycoprotein specifically stored and secreted by neutrophils, which appears to be marker of early differentiation [78, 79]. MCF2L explicitly binds with RAC subfamily of Rho GTPases and may be critical in neutrophilic signal transduction [80, 81]. SPP1 is a major secreted component of pro-inflammatory leukocytes and enhances cellular infiltration and fibrosis of the airways [82, 83, 84]. CFB, increased in both Treated_1 and Treated_2 + cattle when independently compared to Healthy cattle, encodes for the proenzyme of alternative complement. In human sera, CFB has been indicated in early, severe infectious disease, such as that induced by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [85, 86]. Additionally, CFB production is exaggerated in human airway epithelial cells following experimental asthma and rhinovirus challenge [87]. HERC6, IFI6, ISG15, and MX1 encode for proteins induced by type I interferons, such as IFN-α, IFN-β, IFN-ω, and are typically produced by host cells in response to viral infection [88, 89, 90]. Cattle experimentally infected with BRD-associated viruses, such as bovine respiratory syncytial virus (BRSV) and bovine herpesvirus 1 (BHV-1), demonstrate increased expression of these genes during peak clinical presentation, when compared to sham controls [24, 25, 26, 27]. While this appears to be a natural host (bovine) response to viral infection, our previous and current findings, coupled with the work of Sun and colleagues, indicate that anti-viral responses are indicative of early stage and, often, severe naturally-acquired BRD [28, 31]. Previous work indicated strong elevation and co-expression of these genes at arrival in cattle that died of BRD, with predicted interactions with pro-inflammatory cytokine production such as interleukin-6 and tumor necrosis factor-α [31, 91]. Taken together, our findings suggest that cattle that required multiple antimicrobial therapies, and/or died of naturally-acquired BRD, entered the production system with increased anti-viral defense mechanisms without the benefit of activation/elevation of mechanisms for pro-inflammatory regulation or resolution.
To further clarify the relationship between DEGs and disease status, we employed hierarchical clustering and visualization with a gene expression heatmap (Fig. 2). Notably, Healthy and Treated_1 cattle were relatively indistinguishable based on the expression patterns of those 12 DEGs, but Treated_2 + cattle tended to cluster with the association of elevated type I interferon-associated gene expression and decreased SPM and lymphocyte proliferation-associated gene expression (Fig. 2, right side). This finding is similar to results found in our previous multipopulational transcriptome analysis, suggesting that Treated_2 + cattle are more distinct in terms of at-arrival gene expression identification compared to Healthy and Treated_1 cattle [32]. This finding was made more apparent when assessing multilevel modeling and inter-cohort variance, as Treated_2 + cattle tended to be more similar in gene expression compared to Healthy and Treated_1 cohorts, regardless of study-level effect (Figs. 4 and S1). CFB, a gene identified with unequal variance and representing both BRD acquisition and severity, was relatively homogeneous in terms of gene expression within Treated_2 + cattle, but disproportionate within the Healthy and Treated_1 cohorts. Outliers within the Healthy and Treated_1 cohorts tended to shift mean gene expression higher than the median and drove differences in variance observed between the three groups. The observed results of CFB may suggest that some cattle within the Healthy and Treated_1 cohort had subclinical BRD. Further evaluation of variance determined that DSG1, DSP, KRT10, and MGC126945 appeared to be outlier-driven, implying unproportionate differential gene expression across populations, and ruling them out as candidate predictive classifiers for BRD.
Multiclass ROC curve analysis was utilized to generate cutoffs and assess classificational capability for each unique DEG; however, this demonstrated that a single-molecule classification model was not sufficient to capture the dynamic expressional patterns observed across all seven populations (Table 1). Notably, our findings here contrasts with our previous multiclass ROC evaluation to classify disease based on at-arrival differential expression of ALOX15, CFB, MARCO, LOC100335828, and SLC18A2 [32]. While we previously found good to excellent discrimination of Healthy versus Treated_2 + cattle based on the expression levels of single genes in cattle from two independent populations [32], the AUCs resulting from ROC curve analysis in this study were considerably lower when results from all seven populations of cattle were assessed together. For example, in our previous assessment of cattle from two populations, the ROC curve based on differential expression of ALOX15 to distinguish Treated_2 + from Healthy cattle had an AUC of 0.860 [32], while in this study, the AUC was 0.635 (Table 1). The AUCs resulting from ROC curve analysis for each of the seven populations in this study evaluated individually suggests that these molecules may better discriminate BRD in populations with more severe disease, as seen with the DG_17, DG_18, and MH_19 populations compared to VD_17 and PS_19 (Supplemental Tables S1 and S9). For example, in the DG_17, DG_18, and MH_19 populations, the AUCs for the ROC curves differentiating Treated_2 + from Healthy cattle by differential expression of ALOX15 were 0.806, 0.760, and 0.800, respectively; compared to VD_17 and PS_19, for which the AUCs were 0.632 and 0.529, respectively (Supplemental Tables S1 and S9).
Corresponding with our differential expression and hierarchical clustering analyses, Treated_2 + cattle were the most dissimilar and routinely possessed the highest AUC values for genes independently evaluated, when compared to Healthy. Given this, we employed a multivariable ROC curve model for the two most explainable mechanisms (IFN and SPM; Table 1) and developed an expressional decision tree model using the ROC-generated cutoffs from Treated_2 + and Healthy comparisons (Fig. 5). While the combination of ALOX15 and HPGD (SPM) made no discernable difference in performance as compared to evaluation of each gene independently, combining the four IFN genes increased the single-test sensitivity and specificity between Treated_2 + and Healthy cattle. The resulting decision tree possessed relatively low accuracy for discerning Healthy from BRD (51.3%), but possessed excellent accuracy in Treated_2 + categorization (90.0%). While Treated_1 individuals are similar to Healthy in terms of at-arrival gene expression, classificational accuracy (68.5%) is moderately higher than the estimated sensitivity of diagnosing clinical BRD via visual assessment alone (22.0–62.0%) [12, 13]. Furthermore, Treated_1 cattle misclassified as not requiring treatment (Healthy) were sold having higher total average daily gains compared to those classified correctly; this may indicate the ability for this decision tree model to better stratify for disease severity compared to visual assessment and treatment frequency associated with BRD.
This report represents our third evaluation of differential gene expression in whole blood at arrival to predict subsequent BRD treatment in cattle, with progressively larger numbers of cattle from more groups evaluated in each study. Genes that were strongly differentially expressed between Healthy and BRD cattle in our earlier work were not always significantly differentially expressed in subsequent studies. Given the heterogeneity of the cattle within and between these populations, and likely variation in the time of sample collection relative to disease onset, perhaps this finding is no surprise. However, some genes, such as ALOX15 and CFB, have been consistently differentially expressed between Healthy and BRD cattle across our studies, suggesting that further evaluation of these molecules and related pathways may improve prediction and understanding of BRD risk. The influence of the imprecise nature of current BRD assessment systems must also be acknowledged, especially when considering individuals without more definitive clinical endpoints such as mortality. The development of a more accurate method of BRD diagnosis is essential to support efforts to accurately predict the disease.