This study (#6743) was approved by the Institutional Animal Care and Use Committee of Washington State University. A single Idaho dairy provided Dairy Comp 305 (Valley Agricultural Software, Tulare, CA, United States) records and Zoetis CLARIFIDE^® Plus (Zoetis Precision Animal Health, Parsippany, NJ, United States) genotypes of 6,423 Holstein calves. Cattle within the study were born over an eight-year period from a single dairy that milks approximately 2,300 cows. Cattle in the study were from a single herd at a single location in Idaho. Cattle were housed in dry lots with shade and were fed a total mixed ration. All cattle had Dairy Comp 305 records of animal events, health records, and management notes that were used for the study.

The study consisted of a pre-weaned (birth to 60 days of age) calf group and a post-weaned (61–420 days of age) calf group, encompassing the two critical development periods between birth to approximate insemination (14 months). Cases were identified as calves that were treated for BRD or had a recorded respiratory event in the pre-weaned (n = 518) or post-weaned period (n = 2,001). Due to the use of retrospective records, no clinical diagnoses, assessments, or bacteriology/virology were assessed within this study. Diagnosis was solely based upon BRD events, and animals were not removed based upon additional health events, including diarrhea. Control calves were required to remain in the herd for the entirety of the period in which the cases were identified. For the 3,655 healthy, pre-weaned controls, all calves remained and were observed for disease until at least 60 days of age. For the 3,210 healthy post-weaned controls, all calves remained in the herd and were observed for disease for a minimum of 421 days.

Genotypes from Zoetis CLARIFIDE^® Plus (Zoetis Precision Animal Health, Parsippany, NJ, United States) tests were imputed to approximately 620,000 SNPs using a single-step imputation process through Beagle v.4.0 (15) and the ARS-UCD 1.2 assembly (accessed on 8 January 2024).¹ The imputation used a reference population that consisted of roughly 4,800 Holsteins from Washington, Idaho, California, and New Mexico and were genotyped using Illumina BovineHD BeadChip (San Diego, CA, United States). Study animals shared between 29,741 and 53,594 SNPs with the individuals from the reference population and the BovineHD BeadChip. To assess the accuracy of the imputation, previously genotyped animals were also included within the analysis and compared to their known genotypes, determining the final imputation accuracy to be calculated at > 95%. These genotypes were used for the genome-wide association analysis (GWAA) and the gene set enrichment analysis using SNP data (GSEA-SNP).

Prior to GWAA, quality control was completed on the imputed SNPs and the animals genotyped for the pre-weaned and the post-weaned calves. Prior to quality control filtering, the pre-weaned population consisted of 4,171 individuals and 619,410 SNPs and the post-weaned population consisted of 5,211 individuals and 619,410 SNPs. Imputed genotypes were removed if the call rate < 0.9 (n = 6,421 for pre-weaned, n = 6,415 for post-weaned), if they had a minor allele frequency < 0.01 (n = 99,048 for pre-weaned, and n = 98,630 for post-weaned), or if they failed Hardy–Weinberg equilibrium testing with p < 1 × 10⁻¹⁶⁰ (n = 5,985 for pre-weaned, and n = 7,086 for post-weaned). No calves were removed from the pre-weaned or post-weaned group for more than 10% of genotypes failing. After quality control filtering, there were 507,956 SNPs in the pre-weaned heifer analyses and 507,279 SNPs in the post-weaned heifer analyses.

To identify population stratification, a principal component analysis (PCA) was completed for the pre-weaned and post-weaned calves. Clustering by birth year was identified for pre-and post-weaned calves and was included as a covariate in the GWAA. A genomic inflation factor (λ_GC) was calculated to identify the level of population stratification with the use of birth year as a covariate (16). The SNP and Variation Suite (SVS) software version 8.1 (Golden Helix, Bozeman, MT, United States) was used for the GWAA. An efficient mixed-model association eXpedited statistical approach was used to perform the GWAA with an identity-by-state relationship matrix. The general EMMAX statistical approach is defined as y = Χ β + Ζu + ∈ , where y = a n x 1 vector of observed phenotypic values, Χ = an n x f matrix for fixed effects, β = f x 1 vector for the coefficients of fixed effects, Ζ = a matrix containing random effects, u = a vector of random effects with variants of allele substitutions in the population, and ∈ = residual effects (17). Additive, dominant and recessive inheritance models were performed for the GWAA. Associations were established using the Wellcome Trust threshold for uncorrected p-values, with p < 1 × 10⁻⁵ providing evidence for a moderate association and p < 5 × 10⁻⁷ providing evidence for a strong association (18). The proportion of variance explained for each SNP was calculated within SVS for each GWAA. As the proportion of variance explained is not independent for all SNPs, the sum of these effects will exceed 100%. For SNP associated with BRD that were near one another, a locus was characterized when D′ > 0.7 (19, 20). Heritability was estimated within SVS using an AI-REML analysis, which is calculated with a matrix of allele substitution marker effects, otherwise known as a genomic best linear unbiased prediction (GBLUP) (21, 22).

The average haplotype size for the calves in this study was 30,235 bp when calculated using the method of Gabriel et al. (23). Positional candidate genes were identified within one haplotype (± 30 kb) from the associated SNP based on the ARS-UCD 1.2 bovine genome assembly.²

The imputed genotypes were used to conduct the association analysis to identify the SNPs that would serve as representatives for the genes in the gene sets in BioCarta (217 gene sets),³ Gene Ontology or GO (3,147 gene sets),⁴ Kyoto Encyclopedia of Gene and Genomes or KEGG (186 gene sets),⁵ Protein Analysis Through Evolutionary Relationships or PANTHER (165 gene sets),⁶ and Reactome or R (674 gene sets).⁷ The quality control filtered SNPs in the pre-weaned and post-weaned heifer analyses were analyzed and mapped to 21,039 protein-coding genes within the ARS-UCD 1.2 genome assembly. SNPs with the greatest evidence of an association with BRD were used as gene proxies for the GSEA-SNP. Only one SNP was used to represent each gene, noting that a SNP may be a proxy for multiple genes. Genes were ranked by their p-value for their association with BRD from the association analysis. The analysis was run within R v. 3.6.3, using the GenGen v. 1 package (24). An enrichment score for each gene set was computed from the running sum statistics with each gene set receiving a permuted p-value after 10,000 phenotype-based permutations in GenABEL v. 3.6.3 (25, 26). The maximum enrichment score was normalized by the number of genes in the gene set, and those gene sets enriched for BRD had a normalized enrichment score (NES) ≥ 3.0, using a methodology as described by Neupane et al. (27). Leading-edge genes were those that contributed to the peak enrichment score.

There were four loci strongly associated (p < 5 × 10⁻⁷) with BRD in the additive inheritance model, six loci strongly associated with BRD (p < 5 × 10⁻⁷) in the dominant inheritance model, and three loci strongly associated (p < 5 × 10⁻⁷) with BRD in the recessive model in the pre-weaned calves (Table 1 and Figure 1). Two loci, on BTA4 and on BTA7, were strongly associated (p < 5 × 10⁻⁷) with BRD in both the additive and dominant inheritance models. There were 16 additional loci moderately associated (p < 1 × 10⁻⁵) with BRD in the additive model, 11 loci moderately associated (p < 1 × 10⁻⁵) with BRD in the dominant model, and 10 loci moderately associated (p < 1 × 10⁻⁵) with BRD in the recessive model for the pre-weaned calves. The λ_GC values for the pre-weaned BRD population were 0.996 for the additive model, 0.994 for the dominant model, and 0.990 for the recessive model. There were 26 positional candidate genes where an associated SNP fell within an exon or intron of the gene, and an additional 39 positional candidate genes located within a haplotype of a SNP associated with BRD. The estimated heritability for BRD was 0.16 ± 0.02 for the pre-weaned calves.

In the post-weaned calves, 21 loci were strongly associated (p < 5 × 10⁻⁷) with BRD with the additive inheritance model, 18 loci were strongly associated (p < 5 × 10⁻⁷) with BRD with the dominant inheritance model, and 11 loci were strongly associated (p < 5 × 10⁻⁷) with BRD in the recessive inheritance model (Table 2 and Figure 2). Eleven of the loci associated (p < 5 × 10⁻⁷) with BRD in the additive and dominant inheritance models were shared, six loci were associated (p < 5 × 10⁻⁷) in all inheritance models and a single locus associated (p < 5 × 10⁻⁷) with BRD on BTA13 was shared among all three inheritance models. In addition, there were 22 loci moderately associated (p < 1 × 10⁻⁵) with BRD in the additive model, 19 loci associated (p < 1 × 10⁻⁵) with BRD in the dominant inheritance model and 15 loci associated (p < 1 × 10⁻⁵) in the recessive model. The λ_GC values for the post-weaned BRD population were 0.986, 0.995, and 0.996 for the additive, dominant, and recessive models, respectively. There were 56 positional candidate genes where the associated SNP fell within an exon or intron of the gene, and an additional 88 positional candidate genes located within one haplotype of the SNP associated with BRD in the post-weaned calves. The estimated heritability was 0.16 ± 0.02 for BRD in the post-weaned calves.

The Reactome gene set, Phospholipid Metabolism, was enriched (NES = 3.128) for BRD in pre-weaned calves and contained 86 leading edge genes (Table 3). Seven gene sets were enriched (NES ≥ 3) for BRD in the post-weaned calves and consisted of 266 leading edge genes (Table 3). The leading edge gene, p21-activated kinase 1 (PAK1), was present in four of the gene sets enriched for BRD in the post-weaned calves, while FYN proto-oncogene (FYN), cyclin-dependent kinase inhibitor 2B (CDKN2B), cathepsin H (CTSH), casein kinase 2 alpha 1 (CSNK2A1), and integrin subunit beta 1 (ITGB1) were present in three gene sets in post-weaned calves. An additional 90 leading-edge genes were enriched in two gene sets in post-weaned calves. Only a single leading-edge gene, collagen type IV alpha-3-binding protein, also known as ceramide transfer protein or StAR-related lipid transfer protein (COL4A3BP), was shared as a leading-edge gene for BRD in the pre-and post-weaned calves.

This study identified a total of 156 loci across all models for both pre-and post-weaned BRD populations. Previous studies have estimated the heritability of cattle respiratory disease traits to be between 0.04 and 0.28 (14, 28, 29). The estimated heritability within this population was 0.16 for both the pre-and post-weaned calves, which fell within the range of heritabilities previously reported. The estimated heritability of BRD for pre-and post-weaned calves provides encouragement that selection for BRD will be beneficial in reducing disease morbidity and mortality.

A single gene set (phospholipid metabolism) was enriched (NES ≥ 3) for BRD in pre-weaned calves, while seven gene sets were enriched (NES ≥ 3) for BRD in the post-weaned calves (Table 3). None of the gene sets between the two calf groups were shared.

Phospholipid metabolism, enriched for BRD in pre-weaned heifers, is a gene set that includes genes encoding lipids that contain phosphoric acid. Phospholipids are instrumental in maintaining the structure of cell membranes and lipids in foods that provide potential health benefits, such as linoleic acid which is involved with reducing inflammation (92). This pathway contains genes with cellular functions that influence respiratory health (including the production of lipid-based pulmonary surfactants), are involved in cellular structure and in immune signaling during pathogen introduction. The production of mucin, which consists mainly of lipids, is necessary to protect pulmonary epithelium and assist in clearing pathogens (93). Phospholipids present in cell membranes also function to send and receive cellular signals and are key in several lung diseases, including idiopathic pulmonary fibrosis (94). It has been demonstrated that free polyunsaturated fatty acids have an antimicrobial effect by acting on bacterial cell membranes by forming metabolites that affect phagocytosis (95). Whether this mechanism is involved in BRD is unknown.

Neuraminidase 2 (NEU2) is a positional candidate gene and a leading-edge gene within the pre-weaned BRD population in the phospholipid metabolism gene set. NEU2 is one of four members of the sialidase family, which removes sialic acid from glycoconjugates (96). Influenza D is a contributor to BRD (97–99), and influenza D utilizes sialic acid receptors (100, 101). Thus, host NEU2 may function to deplete entry receptors for influenza D virus. Additional work on NEU2 expression in bovine tissues and its influence on viral entry and replication may clarify its relationship to influenza D virus. In addition, NEU2 is upregulated in fibrotic lesions in human and mouse lungs but its role in infectious disease in the lung has yet to be characterized (102, 103).

Of the seven gene sets enriched for BRD in post-weaned heifers, three gene sets (G1 phase, cell proliferation, and regulation of cell proliferation) have functions involving cell proliferation that impacts immune, pulmonary, and epithelial cell production (Table 3). The G1 phase gene set contains nine genes (one leading edge gene) that are important in the first phase of the cell cycle, which also is the phase where external signaling has the greatest ability to pause the cellular proliferation process (104). Viruses, such as influenza A in human host cells, have the ability to arrest the cell cycle at G1 in order to preserve conditions for viral replication (105). Similarly, severe acute respiratory syndrome coronavirus and murine coronavirus can also initiate an arrest between G0 and G1 phases to achieve optimal viral replication (105). Although bacteria do not utilize host cell mechanisms for replication, host immune systems do rely on the proliferation of immune cells to combat bacterial and viral infections. The host’s ability to produce phagocytes and leukocytes is essential in clearing pathogen infections such as those seen in BRD (106).

The 120 genes (6 leading edge genes) involved in the cellular proliferation gene set are essential for basic biological functioning and survival (Table 3). Some of these genes are also involved in the regulation of cell proliferation gene set that contains 95 genes (6 leading edge genes). In humans with tissue damage from emphysema, higher rates of cellular apoptosis of affected cells require higher cellular proliferation to replace the damaged tissue (107). A similar scenario was seen among young rats with chronic respiratory disease where pulmonary cell turnover was twice that of healthy rats of the same age (108). Both the regulation of cellular proliferation through the cell cycle and the proliferation of pulmonary cells themselves highlight research opportunities among cattle to examine the impact of cellular proliferation among cattle lungs with and without BRD.

The remaining enriched gene sets (axon guidance mediated by semaphorins, cysteine type peptidase activity, platelet aggregation and signal transduction by L1) for BRD in post-weaned calves have roles in cell signaling. The gene set axon guidance mediated by semaphorins contains 13 genes (2 leading edge genes) that function in neuronal development and axon formation via semaphorin molecules (109). Beyond axonal guidance, semaphorin proteins function in the immune system, where they moderate T-cell activity (110). Semaphorin proteins regulate neutrophil activation, and migration of immune cells in the inflamed lung as asthma has been associated with unregulated semaphorin expression (111, 112). The semaphorin 3A receptor complex interacts with L1 which is a signal transducing receptor (113).

Signal transduction by L1 is a Reactome gene set that consists of 21 genes (3 leading edge genes). This gene set contains a pathway that directly affects neuronal growth and development. Diseases surrounding L1 largely include neurological conditions such as hydrocephaly and Alzheimer’s disease but is also associated with immune disorders such as fetal and neonatal alloimmune thrombocytopenia, platelet-type bleeding disorder, and erythroleukemia, as well as susceptibility to infectious disease such as dermatitis, anthrax, and west Nile virus (114–116).

The Gene Ontology Cytosine type peptidase activity gene set is a pathway that consists of 12 genes (1 leading edge gene) encoding enzymes that hydrolyze peptide bonds in a polypeptide chain (Table 3). Cysteine peptidase genes are involved in adaptive immune responses by regulating T and B lymphocyte apoptosis (117). Cysteine peptidases facilitate antiviral adaptive immune responses during normal and inflammatory conditions in the lung (118). The dysregulation of peptidases is linked to autoimmune diseases as well as bacterial and viral infections (119). Peptidases can be used by viruses, including coronaviruses, to enter host cells and to aid in viral replication, which has made them potential targets for antiviral treatments (120).

The final gene set is from Gene Ontology and is platelet aggregation that contains only two genes (Table 3). Platelet aggregation has been linked to inflammation, and lung diseases (121). The common function of platelets encompasses blood clotting within internal and external wounds. Moreover, their ability to secrete chemokines recruits immune cells and then binds with the immune cells at locations of pathogen infiltration (122). These platelet-immune cell conglomerates can lead to the phagocytosis of pathogens or further inflammation depending on the cells that bind, and coagulation of platelets also inhibits the spread of bacteria throughout the rest of the body (122). In the lung, platelets serve as the first line of defense to combat alveolar damage from viral and bacterial pathogens in acute respiratory distress syndromes (123). Platelets also have been suggested to have a role in the maintenance of preserving alveolar barriers, both through the production of antioxidant enzymes and restricting alveolar permeability to proteins when in a damaged state (123). Platelets have an important role in the immunological functioning of each individual and also assist in the permeability, inflammation, and maintenance of pulmonary epithelium, which could have important implications in cattle with BRD.

Shared leading-edge genes underscore the important role of those genes in susceptibility to BRD in these calves. There is a single leading-edge gene, collagen type IV alpha-3-binding protein (COL4A3BP), shared between pre-and post-weaned calves. COL4A3BP encodes a ceramide-binding protein, CERT, which is responsible for the transportation of ceramide from its synthesis to its metabolization (124). Ceramides have been linked with cellular membranes and signaling of apoptosis, as they can be converted to ceramide-1-phosphate using ceramide kinase (125). Furthermore, ceramides have been linked to pulmonary inflammation, cystic fibrosis, and emphysema, further highlighting the potential connections between ceramides and BRD (126, 127). A functional connection has been established between ceramides and sphingolipids in the cellular membrane and cystic fibrosis (127). The increase of chronic inflammation within the lungs, influenced via ceramides, increases mucus build-up, inhibits immune responses, and impacts pulmonary structure via cellular apoptosis (127). Pulmonary inflammation and cellular apoptosis are essential elements needed for healthy and properly functioning respiratory systems, and an improper regulation of these mechanisms could lead to an increased risk of contracting diseases such as cystic fibrosis.

In post-weaned calves, there is one leading-edge gene, p21-activated kinase 1 (PAK1) gene on BTA11, that is shared in four of the seven enriched gene sets (NES ≥ 3.0; Table 3). The PAK family of proteins are comprised of six protein members (128). PAK proteins (1–3) play a role in cytoskeleton structure, cell proliferation, and preventing apoptosis (129, 130). In humans, PAK1 proteins have been linked with increased inflammation within the lungs when infected with the COVID-19 virus (131). Whether PAK1 expression also leads to lung inflammation in cattle when exposed to the BRD pathogen, bovine coronavirus, has yet to be established.

Five leading-edge genes, cyclin-dependent kinase inhibitor 2 B (CDKN2B), cathepsin H (CTSH), casein kinase 2 alpha 1 (CSNK2A1), proto-oncogene tyrosine-protein kinase Fyn (FYN), and integrin subunit beta 1 (ITGB1), were shared in three of the seven gene sets enriched (NES ≥ 3.0) for BRD. CDKN2B encodes a protein that inhibits cell cycle progression and has been associated with idiopathic pulmonary fibrosis and amplifies sepsis-induced lung injury (132–134). CTSH has a strong expression in type 2 pneumocytes, and the production of pulmonary surfactant (135). CTSH is differentially expressed in the lungs of individuals at risk for lung adenocarcinoma when macrophages were examined suggesting a role in cellular proliferation in the lung (136). CTSH was also identified as a differentially expressed gene for bovine respiratory disease in a population of Xinjiang calves (137). CSNK2A1 is one of two genes encoding CK2 protein kinases, which function to phosphorylate proteins and can be exploited by viruses who utilize phosphorylated proteins to support viral proliferation (138). CSNK2A1 phosphorylates acid proteins including many transcription factors such as NF-kappa-B, STAT1, CREB1, IRF1, IRF2, ATF1, ATF4, SRF, MAX, JUN, FOS, MYC and MYB as well as proteins involved in immune responses to viral life cycles of Epstein–Barr virus, herpes simplex virus, hepatitis B virus, chronic hepatitis C virus, human immunodeficiency virus, cytomegalovirus and human papillomavirus (139–143). FYN is expressed in T cells and has key roles in the development, selection and maintenance of naïve and peripheral T cells (144, 145). FYN has suggested roles in mice with negative regulation of pulmonary inflammation due to its influence on T-cell signaling (146). Lastly, ITGB1 has been expressed within macrophages and types 1 and 2 pneumocytes, both having implications with pulmonary immunity and maintaining alveolar interface (147). These genes share roles in immune response highlighting opportunities for selection in BRD.

There was also a single gene, SUMO specific peptidase 7 (SENP7), that was identified as a positional candidate and leading-edge genes for BRD in post-weaned calves. The post-translational removal of small ubiquitin-like modifiers (SUMOs) is regulated by the SENP family of genes (148, 149). SENP7 has a role in the formation of heterochromatin during mitosis. Reduced expression of SENP7 can lead to alterations in chromatin structure, as SENP7 is needed for chromatin availability for DNA damage repair (150–152). Genes involved in the regulation of SUMOs have roles in immune cell activation and identifying the presence of, and mounting of, defenses against pathogens (153). Guo et al. (154) found a connection with SENP7 and arthrogryposis multiplex congenita, where one of the symptoms of the fatal disease is early respiratory failure.

Two genes, CERS6 and protein kinase CGMP-dependent 1 (PRKG1), were shared as enriched/associated genes with BRD within both the pre-and post-weaned calves. CERS6, discussed previously due to being highly associated within the post-weaned BRD population, assists in regulating cellular structure and signaling. PRKG1 is located on BTA26 and was a leading-edge gene in pre-weaned heifers and a positional candidate gene in post-weaned heifers. The kinase produced by PRKG1 is responsible for managing smooth muscle relaxation (155). In humans, PRKG1 has functions in bronchodilation and asthma (156–158). Both genes and the role of ceramides and smooth muscle function, provide opportunities for selection for enhanced BRD resistance in cattle.

There was also a shared leading-edge gene, vitamin D receptor (VDR), that was identified within the post-weaned BRD population. VDR was a leading-edge gene within two enriched gene sets (cell proliferation and the regulation of cell proliferation) within the post-weaned BRD population in this study and was also a leading-edge gene among the enriched steroid binding gene set identified by Kiser et al. (27). VDR is responsible for binding vitamin D, which is known to have multiple biological roles including assisting in the regulation of cell cycle control as well as having abilities to repress the expression of T-cells and cytokine producing genes (159, 160). These shared genes highlight the opportunity to identify potential regions for selection for BRD resistance across breed, region, or farm operation.