Comprehensive at-arrival transcriptomic analysis of post-weaned beef cattle uncovers type I interferon and antiviral mechanisms associated with bovine respiratory disease mortality

Host transcriptome response to Mycoplasma bovis and bovine viral diarrhea virus in bovine tissues

BACKGROUND

Mycoplasma bovis is a prominent pathogen associated with respiratory disease in livestock. Respiratory disease in cattle often involves co-infection, where a primary viral infection can weaken the host immune system and thus enhance subsequent bacterial infection. The objective of this study was to investigate changes in the host (cattle) transcriptome during bacterial-viral co-infection. RNA sequencing was done in whole blood cells (WBC), liver, mesenteric lymph node (MLN), tracheal-bronchial lymph node (TBLN), spleen, and thymus collected from Control animals (n = 2), animals infected with M. bovis (MB; n = 3), and animals infected with M. bovis and bovine viral diarrhea virus (BVDV) (Dual; n = 3).

RESULTS

Thymus and spleen had the greatest number of differentially expressed genes (DEGs) out of all tissues analyzed. In spleen, genes involved in maintenance of the extracellular matrix (ECM) including collagen type XV alpha 1 chain (COL15A1), collagen type IV alpha 2 chain (COL4A2), and heparan sulfate proteoglycan 2 (HSPG2) were the most significantly downregulated in Dual compared to Control and MB. In thymus, complement 3 (C3) was a highly significant DEG and upregulated in Dual compared to Control and MB. Interferon alpha inducible protein 6 (IFI6) and interferon-induced transmembrane proteins (IFITM1 and IFITM3), were significantly associated with infection status and upregulated in spleen and thymus of Dual compared to Control and MB.

CONCLUSION

Downregulation of ECM components may cause degradation of the ECM and contribute to increased viral spread due to co-infection. Hyperactivation of complement pathway genes may contribute to damage to the thymus and influence severity of co-infection. Co-expression of IFI6, IFITM1 and IFITM3 across lymphoid tissues may be connected to enhanced pathogenesis in co-infection. These findings suggest co-infection exacerbates disease severity through modulation of ECM components in spleen and complement and coagulation cascades in the thymus. These impacted pathways may underlie thymic atrophy and impaired pathogen clearance due to BVDV and M. bovis co-infection.

Identification of diagnostic biomarkers of and immune cell infiltration analysis in bovine respiratory disease

Background

Bovine respiratory disease (BRD) is a prevalent and costly condition in the cattle industry, impacting long-term productivity, antibioticusage, and global food safety. Thus, identifying reliable biomarkers for BRD is crucial for early diagnosis, effective treatment, and monitoring therapeutic outcomes.

Methods

This study identified differentially expressed genes (DEGs) associated with BRD by analyzing a blood RNA-seq expression dataset associated with BRD, and conducted a Kyoto Encyclopedia of Genes and Genomes (KEGG) approach enrichment and Gene Ontology (GO) annotation analysis on these DEGs. Meanwhile, the key modules related to BRD were screened by weighted gene co-expression network analysis (WGCNA), and the genes in the module were intersected with DEGs. Subequently, least absolute shrinkage and selection operator (LASSO) and random forest (RF) analysis were employed to identify potential biomarkers. Finally, gene set enrichment analysis (GSEA) was performed to explore the potential mechanisms of the identified biomarkers, and their diagnostic significance was assessed using receiver operator characteristic (ROC) curve analysis and real-time fluorescent quantitative PCR (RT-qPCR). In addition, immune cell infiltration in BRD was evaluated using the CIBERSORT algorithm and the correlation between biomarkers and immune cell infiltration was analyzed.

Results

The results showed that a total of 1,097 DEG were screened. GO and KEGG analysis showed that DEGs was mainly enriched in inflammatory response, defense response, Complement and coagulation cascades and Antigen processing and presentation pathways. WGCNA analysis determined that the cyan module had the highest correlation with BRD. A total of 833 overlapping genes were identified through Venn analysis of the differential and WGCNA results. Lasso and RF analyses identified five potential biomarkers for BRD. RT-qPCR testing and data set analysis showed that the expression levels of these five potential biomarkers in nasal mucus and blood of BRD cattle were significantly higher than those of healthy cattle. In addition, ROC curve analysis showed that potential biomarkers had high diagnostic value. GSEA analysis revealed that potential biomarkers are mainly involved in Neutrophil extracellular trap formation, Complement and coagulation cascades, T cell receptor signaling pathway, B cell receptor signaling pathway, Fc gamma R-mediated phagocytosis and IL-17 signaling pathway. The results from the CIBERSORT algorithm demonstrated a significant difference in immune cell composition between the BRD group and the healthy group, indicating that the diagnostic biomarkers were closely associated with immune cells.

Conclusion

This study identified ADGRG3, CDKN1A, CA4, GGT5, and SLC26A8 as potential diagnostic markers for BRD, providing significant insights for the development of new immunotherapy targets and improving disease prevention and treatment strategies.

Comprehensive time-course gene expression evaluation of high-risk beef cattle to establish immunological characteristics associated with undifferentiated bovine respiratory disease

Bovine respiratory disease (BRD) remains the leading infectious disease in beef cattle production systems. Host gene expression upon facility arrival may indicate risk of BRD development and severity. However, a time-course approach would better define how BRD development influences immunological and inflammatory responses after disease occurrences. Here, we evaluated whole blood transcriptomes of high-risk beef cattle at three time points to elucidate BRD-associated host response. Sequenced jugular whole blood mRNA from 36 cattle (2015: n = 9; 2017: n = 27) across three time points (n = 100 samples; days [D]0, D28, and D63) were processed through ARS-UCD1.2 reference-guided assembly (HISAT2/Stringtie2). Samples were categorized into BRD-severity cohorts (Healthy, n = 14; Treated 1, n = 11; Treated 2+, n = 11) via frequency of antimicrobial clinical treatment. Assessment of gene expression patterns over time within each BRD cohort was modeled through an autoregressive hidden Markov model (EBSeq-HMM; posterior probability ≥ 0.5, FDR < 0.01). Mixed-effects negative binomial models (glmmSeq; FDR < 0.05) and edgeR (FDR < 0.10) identified differentially expressed genes between and across cohorts overtime. A total of 2,580, 2,216, and 2,381 genes were dynamically expressed across time in Healthy, Treated 1, and Treated 2+ cattle, respectively. Genes involved in the production of specialized resolving mediators (SPMs) decreased at D28 and then increased by D63 across all three cohorts. Accordingly, SPM production and alternative complement were differentially expressed between Healthy and Treated 2+ at D0, but not statistically different between the three groups by D63. Magnitude, but not directionality, of gene expression related to SPM production, alternative complement, and innate immune response signified Healthy and Treated 2+ cattle. Differences in gene expression at D63 across the three groups were related to oxygen binding and carrier activity, natural killer cell-mediated cytotoxicity, cathelicidin production, and neutrophil degranulation, possibly indicating prolonged airway pathology and inflammation weeks after clinical treatment for BRD. These findings indicate genomic mechanisms indicative of BRD development and severity over time.

Strategies for Bovine Respiratory Disease (BRD) Diagnosis and Prognosis: A Comprehensive Overview

Despite significant advances in vaccination strategies and antibiotic therapy, bovine respiratory disease (BRD) continues to be the leading disease affecting the global cattle industry. The etiology of BRD is complex, often involving multiple microbial agents, which lead to intricate interactions between the host immune system and pathogens during various beef production stages. These interactions present environmental, social, and geographical challenges. Accurate diagnosis is essential for effective disease management. Nevertheless, correct identification of BRD cases remains a daunting challenge for animal health technicians in feedlots. In response to current regulations, there is a growing interest in refining clinical diagnoses of BRD to curb the overuse of antimicrobials. This shift marks a pivotal first step toward establishing a structured diagnostic framework for this disease. This review article provides an update on recent developments and future perspectives in clinical diagnostics and prognostic techniques for BRD, assessing their benefits and limitations. The methods discussed include the evaluation of clinical signs and animal behavior, biomarker analysis, molecular diagnostics, ultrasound imaging, and prognostic modeling. While some techniques show promise as standalone diagnostics, it is likely that a multifaceted approach-leveraging a combination of these methods-will yield the most accurate diagnosis of BRD.

Influence of the At-Arrival Host Transcriptome on Bovine Respiratory Disease Incidence during Backgrounding

Bovine respiratory disease (BRD) remains the leading disease within the U.S. beef cattle industry. Marketing decisions made prior to backgrounding may shift BRD incidence into a different phase of production, and the importance of host gene expression on BRD incidence as it relates to marketing strategy is poorly understood. Our objective was to compare the influence of marketing on host transcriptomes measured on arrival at a backgrounding facility on the subsequent probability of being treated for BRD during a 45-day backgrounding phase. This study, through RNA-Seq analysis of blood samples collected on arrival, evaluated gene expression differences between cattle which experienced a commercial auction setting (AUCTION) versus cattle directly shipped to backgrounding from the cow–calf phase (DIRECT); further analyses were conducted to determine differentially expressed genes (DEGs) between cattle which remained clinically healthy during backgrounding (HEALTHY) versus those that required treatment for clinical BRD within 45 days of arrival (BRD). A profound difference in DEGs (n = 2961) was identified between AUCTION cattle compared to DIRECT cattle, regardless of BRD development; these DEGs encoded for proteins involved in antiviral defense (increased in AUCTION), cell growth regulation (decreased in AUCTION), and inflammatory mediation (decreased in AUCTION). Nine and four DEGs were identified between BRD and HEALTHY cohorts in the AUCTION and DIRECT groups, respectively; DEGs between disease cohorts in the AUCTION group encoded for proteins involved in collagen synthesis and platelet aggregation (increased in HEALTHY). Our work demonstrates the clear influence marketing has on host expression and identified genes and mechanisms which may predict BRD risk.

Differential gene expression in peripheral leukocytes of pre-weaned Holstein heifer calves with respiratory disease

Bovine respiratory disease (BRD) is a leading cause of calf morbidity and mortality, and prevalence remains high despite current management practices. Differential gene expression (DGE) provides detailed insight into individual immune responses and can illuminate enriched pathways and biomarkers that contribute to disease susceptibility and outcomes. The aims of this study were to investigate differences in peripheral leukocyte gene expression in Holstein preweaned heifer calves 1) with and without BRD, and 2) across weeks of age. Calves were enrolled for this short-term longitudinal study on two commercial dairies in Washington State. Calves were assessed every two weeks throughout the pre-weaning period using clinical respiratory scoring (CRS) and thoracic ultrasonography (TUS), and blood samples were collected. Calves were selected that were either healthy (n = 10) or had BRD diagnosed by CRS (n = 7), TUS (n = 6), or both (n = 6) in weeks 5 or 7 of life). Three consecutive time point samples were analyzed for each BRD calf consisting of PRE, ONSET, and POST samples. Nineteen genes of interest were selected based on previous gene expression studies in cattle: ALOX15, BPI, CATHL6, CXCL8, DHX58, GZMB, HPGD, IFNG, IL17D, IL1R2, ISG15, LCN2, LIF, MX1, OAS2, PGLYRP1, S100A8, SELP, and TNF. Comparisons were made between age and disease time point matched BRD and healthy calves as well as between calf weeks of age. No DGE was observed between diseased and healthy calves; however, DGE was observed between calf weeks of age regardless of disease state. Developmental differences in leukocyte gene expression, phenotype, and functionality make pre-weaned calves immunologically distinct from mature cattle, and early life shifts in calf leukocyte populations likely contribute to the age-related gene expression differences we observed. Age overshadows disease impacts to influence gene expression in young calves, and immune development progresses upon a common trajectory regardless of disease during the preweaning period.

Whole blood transcriptome analysis in dairy calves experimentally challenged with bovine herpesvirus 1 (BoHV-1) and comparison to a bovine respiratory syncytial virus (BRSV) challenge

Bovine herpesvirus 1 (BoHV-1), is associated with several clinical syndromes in cattle, among which bovine respiratory disease (BRD) is of particular significance. Despite the importance of the disease, there is a lack of information on the molecular response to infection via experimental challenge with BoHV-1. The objective of this study was to investigate the whole-blood transcriptome of dairy calves experimentally challenged with BoHV-1. A secondary objective was to compare the gene expression results between two separate BRD pathogens using data from a similar challenge study with BRSV. Holstein-Friesian calves (mean age (SD) = 149.2 (23.8) days; mean weight (SD) = 174.6 (21.3) kg) were either administered BoHV-1 inoculate (1 × 10⁷/mL × 8.5 mL) (n = 12) or were mock challenged with sterile phosphate buffered saline (n = 6). Clinical signs were recorded daily from day (d) -1 to d 6 (post-challenge), and whole blood was collected in Tempus RNA tubes on d six post-challenge for RNA-sequencing. There were 488 differentially expressed (DE) genes (p < 0.05, False Discovery rate (FDR) < 0.10, fold change ≥2) between the two treatments. Enriched KEGG pathways (p < 0.05, FDR <0.05); included Influenza A, Cytokine-cytokine receptor interaction and NOD-like receptor signalling. Significant gene ontology terms (p < 0.05, FDR <0.05) included defence response to virus and inflammatory response. Genes that are highly DE in key pathways are potential therapeutic targets for the treatment of BoHV-1 infection. A comparison to data from a similar study with BRSV identified both similarities and differences in the immune response to differing BRD pathogens.

Hematological and gene co-expression network analyses of high-risk beef cattle defines immunological mechanisms and biological complexes involved in bovine respiratory disease and weight gain

Bovine respiratory disease (BRD), the leading disease complex in beef cattle production systems, remains highly elusive regarding diagnostics and disease prediction. Previous research has employed cellular and molecular techniques to describe hematological and gene expression variation that coincides with BRD development. Here, we utilized weighted gene co-expression network analysis (WGCNA) to leverage total gene expression patterns from cattle at arrival and generate hematological and clinical trait associations to describe mechanisms that may predict BRD development. Gene expression counts of previously published RNA-Seq data from 23 cattle (2017; n = 11 Healthy, n = 12 BRD) were used to construct gene co-expression modules and correlation patterns with complete blood count (CBC) and clinical datasets. Modules were further evaluated for cross-populational preservation of expression with RNA-Seq data from 24 cattle in an independent population (2019; n = 12 Healthy, n = 12 BRD). Genes within well-preserved modules were subject to functional enrichment analysis for significant Gene Ontology terms and pathways. Genes which possessed high module membership and association with BRD development, regardless of module preservation (“hub genes”), were utilized for protein-protein physical interaction network and clustering analyses. Five well-preserved modules of co-expressed genes were identified. One module (“steelblue”), involved in alpha-beta T-cell complexes and Th2-type immunity, possessed significant correlation with increased erythrocytes, platelets, and BRD development. One module (“purple”), involved in mitochondrial metabolism and rRNA maturation, possessed significant correlation with increased eosinophils, fecal egg count per gram, and weight gain over time. Fifty-two interacting hub genes, stratified into 11 clusters, may possess transient function involved in BRD development not previously described in literature. This study identifies co-expressed genes and coordinated mechanisms associated with BRD, which necessitates further investigation in BRD-prediction research.

Does swab type matter? Comparing methods for Mannheimia haemolytica recovery and upper respiratory microbiome characterization in feedlot cattle

Background

Bovine respiratory disease (BRD) is caused by interactions among host, environment, and pathogens. One standard method for antemortem pathogen identification in cattle with BRD is deep-guarded nasopharyngeal swabbing, which is challenging, costly, and waste generating. The objective was to compare the ability to recover Mannheimia haemolytica and compare microbial community structure using 29.5 inch (74.9 cm) deep-guarded nasopharyngeal swabs, 16 inch (40.6 cm) unguarded proctology swabs, or 6 inch (15.2 cm) unguarded nasal swabs when characterized using culture, real time-qPCR, and 16S rRNA gene sequencing. Samples for aerobic culture, qPCR, and 16S rRNA gene sequencing were collected from the upper respiratory tract of cattle 2 weeks after feedlot arrival.

Results

There was high concordance of culture and qPCR results for all swab types (results for 77% and 81% of sampled animals completely across all 3 swab types for culture and qPCR respectively). Microbial communities were highly similar among samples collected with different swab types, and differences identified relative to treatment for BRD were also similar. Positive qPCR results for M. haemolytica were highly concordant (81% agreed completely), but samples collected by deep-guarded swabbing had lower amounts of Mh DNA identified (Kruskal–Wallis analysis of variance on ranks, P < 0.05; Dunn-test for pairwise comparison with Benjamini–Hochberg correction, P < 0.05) and lower frequency of positive compared to nasal and proctology swabs (McNemar’s Chi-square test, P < 0.05).

Conclusions

Though differences existed among different types of swabs collected from individual cattle, nasal swabs and proctology swabs offer comparable results to deep-guarded nasopharyngeal swabs when identifying and characterizing M. haemolytica by culture, 16S rRNA gene sequencing, and qPCR.

Supplementary Information

The online version contains supplementary material available at 10.1186/s42523-022-00197-6.

Use of nCounter mRNA profiling to identify at-arrival gene expression patterns for predicting bovine respiratory disease in beef cattle

Background

Transcriptomics has identified at-arrival differentially expressed genes associated with bovine respiratory disease (BRD) development; however, their use as prediction molecules necessitates further evaluation. Therefore, we aimed to selectively analyze and corroborate at-arrival mRNA expression from multiple independent populations of beef cattle. In a nested case-control study, we evaluated the expression of 56 mRNA molecules from at-arrival blood samples of 234 cattle across seven populations via NanoString nCounter gene expression profiling. Analysis of mRNA was performed with nSolver Advanced Analysis software (p < 0.05), comparing cattle groups based on the diagnosis of clinical BRD within 28 days of facility arrival (n = 115 Healthy; n = 119 BRD); BRD was further stratified for severity based on frequency of treatment and/or mortality (Treated_1, n = 89; Treated_2+, n = 30). Gene expression homogeneity of variance, receiver operator characteristic (ROC) curve, and decision tree analyses were performed between severity cohorts.

Results

Increased expression of mRNAs involved in specialized pro-resolving mediator synthesis (ALOX15, HPGD), leukocyte differentiation (LOC100297044, GCSAML, KLF17), and antimicrobial peptide production (CATHL3, GZMB, LTF) were identified in Healthy cattle. BRD cattle possessed increased expression of CFB, and mRNA related to granulocytic processes (DSG1, LRG1, MCF2L) and type-I interferon activity (HERC6, IFI6, ISG15, MX1). Healthy and Treated_1 cattle were similar in terms of gene expression, while Treated_2+ cattle were the most distinct. ROC cutoffs were used to generate an at-arrival treatment decision tree, which classified 90% of Treated_2+ individuals.

Conclusions

Increased expression of complement factor B, pro-inflammatory, and type I interferon-associated mRNA hallmark the at-arrival expression patterns of cattle that develop severe clinical BRD. Here, we corroborate at-arrival mRNA markers identified in previous transcriptome studies and generate a prediction model to be evaluated in future studies. Further research is necessary to evaluate these expression patterns in a prospective manner.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12917-022-03178-8.

Multipopulational transcriptome analysis of post-weaned beef cattle at arrival further validates candidate biomarkers for predicting clinical bovine respiratory disease

Bovine respiratory disease (BRD) remains the leading infectious disease in post-weaned beef cattle. The objective of this investigation was to contrast the at-arrival blood transcriptomes from cattle derived from two distinct populations that developed BRD in the 28 days following arrival versus cattle that did not. Forty-eight blood samples from two populations were selected for mRNA sequencing based on even distribution of development (n = 24) or lack of (n = 24) clinical BRD within 28 days following arrival; cattle which developed BRD were further stratified into BRD severity cohorts based on frequency of antimicrobial treatment: treated once (treated_1) or treated twice or more and/or died (treated_2+). Sequenced reads (~ 50 M/sample, 150 bp paired-end) were aligned to the ARS-UCD1.2 bovine genome assembly. One hundred and thirty-two unique differentially expressed genes (DEGs) were identified between groups stratified by disease severity (healthy, n = 24; treated_1, n = 13; treated_2+, n = 11) with edgeR (FDR ≤ 0.05). Differentially expressed genes in treated_1 relative to both healthy and treated_2+ were predicted to increase neutrophil activation, cellular cornification/keratinization, and antimicrobial peptide production. Differentially expressed genes in treated_2+ relative to both healthy and treated_1 were predicted to increase alternative complement activation, decrease leukocyte activity, and increase nitric oxide production. Receiver operating characteristic (ROC) curves generated from expression data for six DEGs identified in our current and previous studies (MARCO, CFB, MCF2L, ALOX15, LOC100335828 (aka CD200R1), and SLC18A2) demonstrated good-to-excellent (AUC: 0.800–0.899; ≥ 0.900) predictability for classifying disease occurrence and severity. This investigation identifies candidate biomarkers and functional mechanisms in at arrival blood that predicted development and severity of BRD.

Genes and regulatory mechanisms associated with experimentally-induced bovine respiratory disease identified using supervised machine learning methodology

Bovine respiratory disease (BRD) is a multifactorial disease involving complex host immune interactions shaped by pathogenic agents and environmental factors. Advancements in RNA sequencing and associated analytical methods are improving our understanding of host response related to BRD pathophysiology. Supervised machine learning (ML) approaches present one such method for analyzing new and previously published transcriptome data to identify novel disease-associated genes and mechanisms. Our objective was to apply ML models to lung and immunological tissue datasets acquired from previous clinical BRD experiments to identify genes that classify disease with high accuracy. Raw mRNA sequencing reads from 151 bovine datasets (n = 123 BRD, n = 28 control) were downloaded from NCBI-GEO. Quality filtered reads were assembled in a HISAT2/Stringtie2 pipeline. Raw gene counts for ML analysis were normalized, transformed, and analyzed with MLSeq, utilizing six ML models. Cross-validation parameters (fivefold, repeated 10 times) were applied to 70% of the compiled datasets for ML model training and parameter tuning; optimized ML models were tested with the remaining 30%. Downstream analysis of significant genes identified by the top ML models, based on classification accuracy for each etiological association, was performed within WebGestalt and Reactome (FDR ≤ 0.05). Nearest shrunken centroid and Poisson linear discriminant analysis with power transformation models identified 154 and 195 significant genes for IBR and BRSV, respectively; from these genes, the two ML models discriminated IBR and BRSV with 100% accuracy compared to sham controls. Significant genes classified by the top ML models in IBR (154) and BRSV (195), but not BVDV (74), were related to type I interferon production and IL-8 secretion, specifically in lymphoid tissue and not homogenized lung tissue. Genes identified in Mannheimia haemolytica infections (97) were involved in activating classical and alternative pathways of complement. Novel findings, including expression of genes related to reduced mitochondrial oxygenation and ATP synthesis in consolidated lung tissue, were discovered. Genes identified in each analysis represent distinct genomic events relevant to understanding and predicting clinical BRD. Our analysis demonstrates the utility of ML with published datasets for discovering functional information to support the prediction and understanding of clinical BRD.

Integrated Network Analysis to Identify Key Modules and Potential Hub Genes Involved in Bovine Respiratory Disease: A Systems Biology Approach

Background: Bovine respiratory disease (BRD) is the most common disease in the beef and dairy cattle industry. BRD is a multifactorial disease resulting from the interaction between environmental stressors and infectious agents. However, the molecular mechanisms underlying BRD are not fully understood yet. Therefore, this study aimed to use a systems biology approach to systematically evaluate this disorder to better understand the molecular mechanisms responsible for BRD. Methods: Previously published RNA-seq data from whole blood of 18 healthy and 25 BRD samples were downloaded from the Gene Expression Omnibus (GEO) and then analyzed. Next, two distinct methods of weighted gene coexpression network analysis (WGCNA), i.e., module-trait relationships (MTRs) and module preservation (MP) analysis were used to identify significant highly correlated modules with clinical traits of BRD and non-preserved modules between healthy and BRD samples, respectively. After identifying respective modules by the two mentioned methods of WGCNA, functional enrichment analysis was performed to extract the modules that are biologically related to BRD. Gene coexpression networks based on the hub genes from the candidate modules were then integrated with protein-protein interaction (PPI) networks to identify hub-hub genes and potential transcription factors (TFs). Results: Four significant highly correlated modules with clinical traits of BRD as well as 29 non-preserved modules were identified by MTRs and MP methods, respectively. Among them, two significant highly correlated modules (identified by MTRs) and six nonpreserved modules (identified by MP) were biologically associated with immune response, pulmonary inflammation, and pathogenesis of BRD. After aggregation of gene coexpression networks based on the hub genes with PPI networks, a total of 307 hub-hub genes were identified in the eight candidate modules. Interestingly, most of these hub-hub genes were reported to play an important role in the immune response and BRD pathogenesis. Among the eight candidate modules, the turquoise (identified by MTRs) and purple (identified by MP) modules were highly biologically enriched in BRD. Moreover, STAT1, STAT2, STAT3, IRF7, and IRF9 TFs were suggested to play an important role in the immune system during BRD by regulating the coexpressed genes of these modules. Additionally, a gene set containing several hub-hub genes was identified in the eight candidate modules, such as TLR2, TLR4, IL10, SOCS3, GZMB, ANXA1, ANXA5, PTEN, SGK1, IFI6, ISG15, MX1, MX2, OAS2, IFIH1, DDX58, DHX58, RSAD2, IFI44, IFI44L, EIF2AK2, ISG20, IFIT5, IFITM3, OAS1Y, HERC5, and PRF1, which are potentially critical during infection with agents of bovine respiratory disease complex (BRDC). Conclusion: This study not only helps us to better understand the molecular mechanisms responsible for BRD but also suggested eight candidate modules along with several promising hub-hub genes as diagnosis biomarkers and therapeutic targets for BRD.