In recent years, bovine RNA-Seq data has been utilized to identify biomarkers that are predictive of or associated with BRD. Previous studies have utilized tissue samples, such as lung, lymph node, and tonsil, at the time of peak clinical infection following challenge with individual BRD pathogens13,14,15. These studies provided important foundational efforts for identifying regulatory networks and host-pathogen interactions related to immunological defense at the time of peak clinical BRD. However, pathogen challenge models are limited in the degree to which they represent naturally occurring BRD, because they do not recapitulate the complex multifactorial nature of the disease47,48,49. In contrast, gene expression patterns have been analyzed in arrival blood from post-weaned beef cattle that later developed BRD16,17,18. These investigations of naturally occurring BRD assure the biological context that is necessary for gene expression analyses to identify altered biological systems for clinical BRD prediction. By testing large populations of BRD susceptible cattle, this approach allows differentiation of pervasive changes in biological systems that segregate with BRD prediction and staging, from other dynamic environmental, and host factors that also alter biological systems but do not segregate with BRD. Using this approach, altered immunological and metabolic pathways, including those important to host-pathogen interactions, are recognized by their shared patterns in cattle that develop BRD. Patterns that influence clinical BRD outcomes and support the discovery of novel pathophysiology, candidate biomarkers and therapeutic modalities can also be realized 50,51.
The identification of differentially expressed genes and enriched genomic mechanisms within whole blood samples supports the discovery and use of candidate biomarkers and novel therapeutic modalities. Whole blood is an easily obtainable, non-invasive sample type that represents transient biological occurrences at distinct physiological sites52,53. By analyzing transcriptomes in blood collected at arrival from two distinct populations of cattle, this investigation 1) corroborated our prior findings that at-arrival expression profiles of MARCO, CFB, MCF2L, ALOX15, LOC100335828 (CD200R1), and SLC18A2 were common to cattle from two distinct populations that went on to develop BRD, 2) demonstrated using ROC curves that the at-arrival expression of these genes have predictive potential in classifying BRD acquisition and severity, and 3) derived novel genomic information relevant to BRD severity, related to leukocyte activity and airway epithelium.
A limitation of this study is the many factors that varied between individual cattle and between the populations which contribute to biological variation in the transcriptomes of individual cattle. In this regard, our principal component analysis identified population (Year) in PC3 as a significant and identifiable source of variation in gene expression (Fig. 3B). Our populations lacked information regarding treatment/vaccination of individual animals prior to arrival. Additionally, all individuals in this study were of varying breeds, were likely to be at differing time points in the disease spectrum of BRD, and were likely to have a host of other factors that varied between individuals prior to their arrival. Relevant to this variation is our overarching objective, which is to identify arrival gene expression that predicts BRD outcome in commercial beef cattle within 28 days following arrival. To achieve this objective, it is essential that the very factors causing variation in the transcriptomes of beef production systems, both recognized and implicit, are included in the experimental model. Only within this biological context is it possible for experimental techniques, such as PCA that was employed within this investigation, to be used to identify correlations between gene expression, disease, and both known and implicit sources of variation. An additional relevant limitation of this study is that it is underpowered for its overarching goal of identifying genes whose differential expression can be universally applied to all at arrival beef cattle to predict BRD within 28 days. However, our approach is an early, foundational event in support of this goal. Our demonstrated ability from ROC curves analysis, to predict calves at arrival that will develop BRD within 28 days, provides proof of concept for the eventual validation of predictive at-arrival biomarkers for clinical BRD. These findings also highlight the necessity to expand our approach to include cattle populations that will assure the appropriate power (i.e., number of animals and of unique herds) to characterize gene expression that predicts developing BRD in the face of gene expression variations that are characteristic of these populations.
Several studies have shown that production and economic loss increases with an increase in frequency of treatment and earlier timing of initial BRD treatment11,54,55,56. To align to this insight, and predict differences in future BRD severity, we stratified the severity of the BRD phenotype in our experiment. Cattle were separated into BRD severity cohorts based upon frequency of antimicrobial treatment, clinical assessment scores, and BRD-associated mortality. We identified 132 unique DEGs between the three disease cohorts. The increased DEGs identified in treated_1 cattle, when compared to both healthy and treated_2+, were largely involved with three major innate functions: neutrophil recruitment and degranulation, antimicrobial peptide production, and cellular cornification/keratinization. Increases in neutrophil recruitment/degranulation and antimicrobial peptide production traced primarily to increased differential expression of CATHL1/2/3, CAMP, DSP, PRG3, AZU1, CTSG, CD177, MPO, LTF, and NGP. These products are important for host antimicrobial defense, particularly involving leukocytic interactions within the airways57,58,59,60, and several of these gene products have been directly identified in cattle that were experimentally infected with BRD agents13,14,57.
Cattle produce all known mammalian classes of cathelicidins, compared to humans which only produce one, and research has primarily focused with their effects on bacterial pathogens61,62,63. Antimicrobial peptides, particularly cathelicidins, as well as CD177-mediated neutrophilic response in cattle, have been shown to be effective at mitigating Gram-negative bacterial infections and may be associated with modulating the cytotoxic/apoptotic host responses through the decrease of extracellular IL1A availability58,61,64. Dysregulation of these cytotoxic/apoptotic pathways is linked with poor clinical outcomes57,65,66. The increased cellular cornification/keratinization functions in treated_1 cattle versus healthy and treated_2+, identified through GO and pathway enrichment analyses, were impacted by increased expression of KRT5/10/14/25, DSG1, DSP, CDSN, and SPINK5 in this group. While keratinization is primarily considered a modification of epidermal skin cells, this process has been shown to be an anti-apoptotic process of host barrier defense and structural repair in the lower airways related to infectious respiratory disease67,68.
To reduce dimensionality of the data and enable correlations between known sources of variability in our experiment and variability in gene expression to be identified, we performed principal component analysis that included both gene expression data and metadata for known variables for each animal. Component loadings identified positive correlation to disease severity and negative correlation to average daily gain in PC1, which were of moderate effect and statistically significant (Fig. 3B; FDR <0.05). Accordingly, PC1, which by definition captures the greatest variability in the data, primarily measures the degree of weight loss and disease severity in BRD-affected cattle. This aligns to the known correlation between weight loss over time and disease in beef cattle54,56. In addition, CFB, a gene with the greatest positive influence on PC1 variation and a positive correlation to measures of disease severity using PCA, was also identified as acceptable discriminator of cattle that developed more severe disease (severe_2+) from cattle that remained healthy, based upon ROC analysis (Fig. 4; AUC= 0.769).
CFB, which was increased in both treated_1 and treated_2+ cattle when compared to healthy (Supplemental Table S4), encodes for complement factor B, the major acute phase protein needed for activation of the alternative pathway of complement. CFB is a pro-inflammatory molecule, produced by type II alveolar epithelial cells69,70. CFB can directly stimulate monocyte and B lymphocyte response to viral and bacterial respiratory pathogens69,70,71. CFB has been consistently identified as significantly expressed in BRD-afflicted cattle, both from pathogen challenge models and at-arrival sampling in naturally occurring BRD13,14,15,16,17,18. CFB is secreted by and activates M1 macrophages72,73,74. Our network analysis (Fig. 5; blue cluster) and IPA (Table 1 & 2) indicates that treated_2+ cattle, compared to treated_1 cattle, possess increased pro-inflammatory activity, which is driven by IL1A and FADD. This is made evident by the increased expression in macrophage-specific receptors MARCO and CD163, and decrease in the inhibitor glycoprotein CD200R1, when compared to both healthy and treated_1 cattle. Evidence suggests that MARCO, CD163, and CD200R1 expression, induced by ligand binding to these receptors, causes a pro-inflammatory response by tissue macrophages, particularly in the lung, and may lead to an enhanced IL-1/IL-6 response75,76,77.
Through our previous and current research, we identified the differential expression of several genes involved in lipid transport and anti-inflammatory modalities16. Several genes that correspond to multidrug resistance-associated protein 4 (ABCC4/MRP4) were significantly decreased in treated_2+ cattle compared to both healthy and treated_1. MRP4 is a transporter protein found in multiple cell types and is directly involved in lipid molecule transport, such as prostaglandins and eicosanoids78,79. When prostaglandin E2 (PGE2) is produced via cyclooxygenase-2 (COX-2), MRP4 exports PGE2 to the extracellular space for pro-inflammatory signaling78,80,81. However, MRP4 also mediates the efflux of eicosanoids and its expression is coupled with an increase in ALOX15 expression. ALOX15 encodes for the lipid peroxidizing enzyme arachidonate 15-lipoxygenase, expressed by airway epithelium, circulating reticulocytes, macrophages, mast cells, and eosinophils82,83,84. ALOX15 is required for the biosynthesis of specialized proresolving mediators (SPMs), such as resolvins and lipoxins83,85. These SPMs are important anti-inflammatory mediators derived from arachidonic acid and polyunsaturated fatty acids, and are responsible for metabolizing and suppressing the effects of prolonged inflammatory mediator signaling, including PGE278,83,85. ALOX15 production is implicated as a mitigating factor in a wide variety of inflammatory diseases in humans86,87. Regarding BRD, ALOX15 has been found to be differentially expressed in animals challenged with single pathogens; importantly, those animals were all determined to exhibit mild clinical disease at time of sampling13,14,15. Additionally, at-arrival blood transcriptomes from this and previous research have shown cattle that remain healthy exhibit increased expression of ALOX1516. Thus, it can be hypothesized that cattle in these production settings are undergoing some form of cellular stress and subsequent pro-inflammatory signaling, but ALOX15 upregulation and SPM production may be protective against severe clinical BRD.
We selectively evaluated the performance of six candidate mRNA biomarkers and their ability to predictively classify severity-based cohorts with ROC curves and calculated AUCs (Fig. 4). The genes evaluated in this analysis were identified as differentially expressed in both this current and our previous studies (MARCO, CFB, MCF2L, ALOX15, LOC100335828 (CD200R1), and SLC18A2)16. This approach demonstrated excellent discernment, at arrival and prior to the onset of clinical signs, or cattle that would develop severe BRD (treated_2+) within 28 days of arrival. Accordingly, these genes have inherent prognostic value for accurately identifying animals that require the highest frequency of treatment, yield the lowest weight gains overtime, and/or succumb to BRD-associated mortality. The ability to accurately predict cattle that are most likely to develop severe BRD allows for precise management and antimicrobial treatment protocols, reducing the need to medicate beef cattle that are less likely to develop BRD during their production phases. Furthermore, predicting cattle that require multiple antimicrobial treatments, or eventually succumb to BRD, has significant economic impact. Several studies have established that beef cattle that receive multiple BRD treatments over their production phases yield lower carcass grades and weights, increased management and treatment cost per day, and overall reduced economic returns overtime88,89,90. Future research is necessary for evaluating these genes as prognostic indicators in larger and more complex beef cattle populations.