Association of ARRDC3 and NFIA variants with bovine congestive heart failure in feedlot cattle

Ethical statement

The experimental design and procedures used during this research project were reviewed and approved by the Institutional Animal Care and Use Committee of the University of Nebraska-Lincoln as previously described (Experimental Outline numbers 139 and 1172 [ Workman et al., 2016; Heaton et al. 2019]). The animals were privately owned by commercial feeding operations, and the owners and management approved the use of animals for this study. In every instance, all efforts were undertaken to reduce animal suffering. In addition, animal welfare for these facilities was in accordance with the National Cattlemen’s Beef Association’s Beef Quality Assurance Feedyard Welfare Assessment program (http://www.bqa.org/).

Animals and study design

Paired samples from 102 affected feedlot cattle and their 102 unaffected matched pen mate controls were identified from approximately 140,000 total animals collected from four feedlots in Nebraska and Wyoming at elevations ranging from 1,163 to 1,280 m, in a 16-month period spanning January 2017 to April 2018 (Heaton et al., 2019). Briefly, end-stage BCHF cases were identified and euthanized by feedlot personnel based on clinical presentation as previously described. The case definition included multiple clinical signs specific and non-specific to BCHF as described previously in detail (Heaton et al., 2019). Cases were enrolled in the study only if there was a postmortem presumptive diagnosis of congestive heart failure at necropsy. Unaffected pen mate controls were matched for source, arrival date, gender, and coat color. The 102 case-control pairs originated from 30 different sources with the largest single source contributing 32 pairs (Heaton et al., 2019). The pen sizes ranged from 100 to 300 animals, and control animals were selected based on their willingness to be moved to the sampling area by riders on horseback.

None of the unaffected matched pen mate animals developed clinical BCHF signs prior to their processing into beef. Overall, the group of 204 calves was 93% solid black, 70% castrated males, and none had visible horns (Heaton et al., 2019). Most of the 204 calves were from known, well-managed herds that use American Angus genetics. However, breed information was not available for all animals. Tissue samples for genomic DNA isolation included V-shaped ear notches and EDTA whole blood. The ear tissue was desiccated with granular NaCl in the field and stored at -20°C upon return. The plasma and cellular fractions of EDTA whole blood were separated and frozen in liquid nitrogen en route and stored at -80°C upon return. Postmortem whole hearts used for gross morphological comparisons representing unaffected cattle were collected during federally-inspected beef processing at the USMARC abattoir from 21 purebred American Angus heifers raised and fattened at an elevation of 578 m.

A panel of 96 unrelated individuals from popular U.S. beef breeds (USMARC Beef Diversity Panel version 2.9 [MBCDPv2.9]) was used for genetic comparisons with case-control animals (Heaton et al., 2016). The panel design was based on a previous set of commercially-available sires from 19 breeds with minimal pedigree relationships (Heaton et al., 2001). Pedigrees were obtained from leading suppliers of U.S. beef cattle semen and analyzed to identify the least related individuals for inclusion. On the basis of the number of registered progeny, the breeds were estimated to represent more than 99% of the germplasm used in the U.S. beef cattle industry. Based on pedigrees, the panel members contained more than 187 unshared haploid genomes (of 192 possible) and allowed a 95% probability of detecting any allele with a frequency greater than 0.016.

Genotype scoring, SNP filtering, and population substructure evaluation

Unless otherwise indicated, reagents were molecular-biology grade. DNA from ear notches was extracted by standard phenol/chloroform procedures, dissolved in a solution of 10 mM TrisCl, 1 mM EDTA (TE, pH 8.0) and stored at 4°C. (Heaton et al., 2008). SNP genotypes for variants on the BovineHD BeadChip array were scored at GeneSeek Inc. (Lincoln, NE, USA), according to manufacturer’s instructions (Illumina, Inc., San Diego, CA, USA). The unfiltered map and ped files from the 102 BCHF case-control animals (Files S1 and S2, Extended data) were analyzed for association with BCHF with PLINK v1.90b6.5 64-bit (13 September 2018) software (Chang et al., 2015). The pairs (1 to 102) were designated in the PLINK family ID field. SNPs were removed for having a minor allele frequency less than 0.05 (--maf option). This cutoff was conservative since the frequency of the affected individuals was set at 0.50 a priori in this matched case-control design. The missing genotype rate threshold was set at 0.05 (--geno option) and the Hardy-Weinberg equilibrium test with the mid-p adjustment threshold was set to 0.0001 (--hwe midp option) (Graffelman and Moreno, 2013). These filters were applied for 204 samples of which 142 were listed as males (steers) and 62 were listed as females (heifers). The animal genders were tested with the --check-sex option and a control animal was identified as being misclassified as a male, leaving 141 steers and 63 heifers, and pair 21 mismatched for sex with the case being a steer and the control being a heifer (Files S3 and S4, Extended data).

Since unrecognized population stratification among individuals may obscure association signals and lead to incorrect findings (Deng, 2001; Cardon and Palmer, 2003), multidimensional scaling (MDS) was used with a matrix of pairwise identity-by-state (IBS) distance to visually evaluate the similarity between individuals and pairs. Using PLINK software, unfiltered input map and ped files from the 102 BCHF case-control animals (Files S1 and S2, Extended data) were merged with likewise unfiltered input map and ped files from the USMARC Beef Cattle Diversity Panel version 2.9 (MBCDP v2.9) genotyped with BovineHD BeadChip array (Heaton et al., 2016) (Files S5 and S6, Extended data). Once merged, the data was filtered as described above to produce a new map and ped files (Files S7 and S8, Extended data). MDS data was produced with PLINK software v1.9 and plotted with Microsoft Office 365 Excel software. The IBS/IBD computation was run on the filtered SNP sets with the --genome option. The computation excluded 12,591 variants on non-autosomes from the IBD calculation. Plotting the --cluster output values in the first and second dimensions (i.e., C1 values against C2) gave a scatter plot in which each point is an individual; the two axes correspond to a reduced representation of the data in two dimensions, which facilitates cluster identification. Standard classical (metric) multidimensional scaling was used.

McNemar’s test with SNP genotypes

The McNemar’s test (McNemar, 1947) is a statistical test for paired nominal data. Unlike a traditional case-control GWAS where the scope of inference is over two large groups of individuals, a McNemar’s test analysis for a case-control GWAS uses many pairs of individual animals. In addition, McNemar’s test analysis is not formally available as a convenient option in commonly used, open-source GWAS toolsets that can analyze 560,000 SNPs such as PLINK (v1.90b6.5). The Cochran-Mantel-Haenszel (CMH) test (Mantel and Haenszel, 1959) is available in PLINK and can be run as a specialized case with two clusters; however, it has limited options and access to essential McNemar raw data is not available. Thus, a custom workflow was developed with a programming and numeric computing platform (MatLab v.2021.2, MathWorks, Inc., Natick, MA) to produce the critical tables of binomial outputs and McNemar statistics for the filtered SNP genotypes.

In McNemar’s test, informative pairs are those which contain exactly one animal with the hypothetical risk allele(s), i.e., Table 1, quadrants b and c. The ratio of these discordant pairs (b/c) provides evidence to reject the null hypothesis. At each biallelic SNP locus, alleles were tested separately in three scenarios as a potential disease risk factor: 1-copy, 1- or 2-copies, and 2-copies, producing six sets of McNemar’s test statistics. The presence of a potential risk allele(s) was scored for both animals in every matched pair and used to populate the McNemar’s 2 × 2 contingency table for a given SNP. The presence (+) or absence (-) of a possible genetic risk factor was assigned to each animal and the matched pair was then classified by one of four possible binomial outcomes (Table 1, quadrants). A genetic “risk” allele was indicated in the 102 matched pairs when quadrant b was greater than quadrant c (i.e., odds ratio (OR) > 0). Conversely, a “protective” allele was indicated when quadrant c was greater than quadrant b (i.e., OR < 0).

McNemar’s test analysis for a GWAS imposes statistical constraints at the SNP level that are not present in conventional GWAS. The McNemar’s test p-values for rejection of the null hypothesis were computed with the binomial distribution. This is unlike a conventional retrospective case-control GWAS, where the theoretical distribution of p-values may be modeled with Fisher’s exact test statistics (or derivative). The binomial distribution generated discrete p-values (e.g., exact p and mid p) at each SNP site that were dependent on the discordant pairs observed (i.e., b + c and b/c) for each inheritance model (e.g., 1-copy, 1- or 2-copies, and 2-copies). Thus, p-values were computed de novo at each locus. To the best of our knowledge, the theoretical genome-wide distribution of McNemar’s p-values at the SNP level has not been previously described in the GWAS context. Briefly, input files were PLINK map and ped files with the ped file modified to include the pair identifier in the first column, and the rows sorted in descending pair order from 1 to 102 with the unaffected animal listed first in the pair. The case and control phenotype were coded in the sixth column, with the unaffected control animal coded with 1 (PHENO = 1) and the affected case animal with a 2 (PHENO = 2). The analysis evaluated each SNP for possible risk alleles in the context of McNemar’s test as described in README.md file in the GitHub code repository McNemarsSNPAnalysis. Alternatively, the analysis pipeline can be run online within a copy of the McNemar’s SNP Analysis Compute Capsule at CodeOcean. The PLINK map and ped files contained 141 steers and 63 heifers, and SNPs assigned to chromosome 0 through 32 (i.e., included unassigned SNPs, chromosome X, Y, and mitochondrial DNA; Files S9 and S10, Extended data). In the first stage of the workflow, biallelic SNP heterozygous animal genotypes were converted to a one-letter International Union of Pure and Applied Chemistry and International Union of Biochemistry and Molecular Biology (IUPAC/IUBMB) ambiguity codes used for nucleotides (i.e., R = a/g, Y = c/t, M = a/c, K = g/t, S = c/g, W = a/t (‘Nomenclature Committee for the International Union of Biochemistry [NC-IUB]. Nomenclature for incompletely specified bases in nucleic acid sequences. Recommendations 1984’, 1986) and homozygous animal genotypes were converted to a one-letter code of the same letter (i.e., A, C, T, G).

The McNemar mid-p test was used, since it is an improvement over the exact conditional test (Fagerland et al., 2013). Briefly, the McNemar mid-p-value was calculated with equation 1:

The quantile-quantile (Q-Q) plots of the distribution of the test statistics are challenging to produce because McNemar’s test data are discrete for a given number of total pairs and the fraction of discordant pairs ((b + c)/n, Table 1). This requires the computed theoretical expected distribution of p-values to mirror the dependency of the observed distribution of p-values on the fraction of discordant pairs. Thus, for a given informative pair (b + c) value, theoretical p-values in Q-Q plots were computed by multiplying the number of observed p-values by the predicted normalized binomial distribution of all possible b and c values, as defined by the p-value metric used, e.g. exact p or mid p. The fraction of discordant pairs is also influenced by minor allele frequency, as the opportunity for the observation of discordant pairs diminishes with frequency.

Systematic use of q-values in genome-wide tests provides a balance between the number of true and false positives while being automatically calibrated and readily interpreted (Storey and Tibshirani, 2003). The q-value of a genome-wide data set is the expected proportion of false positives incurred (Storey and Tibshirani, 2003). For instance, a q-value threshold of 5% results in the false discovery rate (FDR) of 5% of the significant features being truly null. The mafdr function in the Matlab Bioinformatics Toolbox (version 2021.2) was used to compute the q-values from input p-values according to procedures previously described (Storey, 2002). The tuning parameter, λ, was used to estimate the a priori probability that the null hypothesis is true (equation 2):

The six sets of McNemar’s pair classifications and test statistics were tallied for each SNP allele and exported to a text file for importation into spreadsheet software (Excel Version 1902, Build 11328.20420). A spreadsheet environment allowed for the SNP set to be used in calculations for McNemar’s statistics such as odds ratios (OR), p-values, and McNemar Chi-squared statistics. After McNemar’s testing, SNPs were removed due to their assignment to chromosomes 0, 31, and 32. These SNPs represent markers on unmapped chromosomes, SNPs on chromosome Y, and mitochondrial SNPs. SNPs were also removed if there were no informative pairs resulting from their pairwise analysis (i.e., McNemar quadrants b + c = 0).

The upper and lower CI₉₅ for the OR were calculated in equations 3 and 4:

The spreadsheet format allowed for custom labeling of column headers, row sorting, and color-coded heat mapping. For example, it was useful to sort each of the six sets of data by their best McNemar’s p-value in descending order and color code the rows in blocks of 100 (e.g. purple, red, yellow, dark green, light green). This helped visually identify the top 500 best markers in each data set. When sorted by chromosome and position, this produces an obvious visual symmetry due to the biallelic nature of SNPs and the dichotomous phenotype.

The spreadsheet format also allowed for manual production of Manhattan plots from McNemar’s data for three inheritance models of the two alleles for each SNP. Briefly, the ARS-UCD1.2 reference assembly SNP coordinates from the PLINK map file (i.e., chr and position) were converted to a linear genome position for the x-axis of the Manhattan plot so that every SNP had a new relative coordinate on a whole genome scale from 1 to 2,651,083,094 bp. This was accomplished by adding the position value of the most terminal SNP on chromosome 1 to each of the SNP position values on chromosome two. This was repeated for each consecutive chromosome (1 to 29 and X). The Bonferroni correction for multiple testing was used at the 0.05 significance level (α) as a conservative genome-wide estimate of significance (Johnson et al., 2010) (Equation 5):

Power analysis

Sample size power analysis for the McNemar’s test was performed with PASS 2019 Power Analysis and Sample Size Software, version 19.0.2 (NCSS, LLC., Kaysville, Utah, USA). The sample size of 102 matched pairs achieved 95% power to detect an OR of 8.51 and 2.94 with a two-sided McNemar’s test and a corrected significance level of 0.01 (α) within the range of 25% and 70% discordant pairs, respectively. Power is the probability of rejecting a false null hypothesis, while α is the probability of rejecting a true null hypothesis.

CMH test

McNemar’s test analysis was compared to that obtained with the CMH association statistic (Mantel and Haenszel, 1959) in PLINK v1.90b6.5. The 102 matched pairing assignments were coded in the fam file in the third column to identify 102 “clusters” (e.g., Pair₁, Pair₂, Pair₃, …Pair₁₀₂) for use with the CMH association statistic. The --adjust option was also used to produce a sorted list of CMH association results. The Q-Q plot of the distribution of the test statistics was produced with the qqman R package on CRAN: (Turner, SD. 2014. bioRxiv) with 64-bit R version 3.5.1 (2018-07-02).

Testing BCHF cases and controls for bovine viral diarrhea virus (BVDV) infection

Blood samples from BCHF cases and controls were tested for infection with BVDV, a possible confounding disease condition of feedlot cattle that is easily detected. However, since these animals were not serially sampled, a positive test result would not distinguish between acute and persistent infection. Blood samples for the first 13 case-control pairs were collected as EDTA whole blood, stored at -80^oC, and used for BVDV testing. Blood samples for the remaining 89 case-control pairs were processed by first separating the plasma by centrifugation, and storing it separately at -80^oC prior to its use in BVDV testing. EDTA blood or plasma was pooled into groups of five or 10, respectively, and total RNA was extracted with a mono-phasic solution of phenol and guanidinium thiocyanate (TRIzol LS, Invitrogen-Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s specifications. BVDV genomic RNA was detected with a previously published, real-time assay that used quantitative reverse transcription polymerase chain reaction (RT-qPCR) primers (5′-GGGNAGTCGTCARTGGTTCG-3′; 5′-GTGCCATGTACAGCAGAGWTTTT-3′) and probe (5′-6-FAM-CCAYGTGGACGAGGGCAYGC-BHQ1-3′) (Mahlum et al., 2002). The assay reagents were purchased as a commercial kit (Qiagen OneStep RT-PCR, Germantown, MD) and used according to the manufacturer’s instructions. Briefly, cyclic amplification was conducted in a 25 μL reaction containing 5 μL of extracted RNA, 4.5 mM magnesium chloride, 400 μM dNTP, 0.4 μM of each forward and reverse primer, 0.2 μM probe, and 1 μL of an enzyme mix containing reverse tran-scriptase and a “hot-start” Taq polymerase. Cycling conditions included reverse transcription at 50°C for 30 minutes, inactivation of the reverse transcriptase enzyme and activation of the Taq polymerase at 95°C for 15 minutes, and 40 cycles of 94°C for 30 seconds, 55°C for 60 seconds, and 72°C for 60 seconds. Positive, negative, no-template, and extraction control samples were included in each group of samples tested.

SNP genotype assays for ARRDC3 and NFIA

Custom SNP genotype tests were developed on a commercially available genotyping platform that uses “padlock” oligonucleotide probes, combined with fluorescently-labeled probes and isothermal rolling circle amplification technology, for the qualitative detection of alleles (MatMaCorp, Lincoln, NE, USA). Padlock probes are long oligonucleotides whose ends are complementary to adjacent target sequences. The primary target SNPs included ARRDC3 C182Y (ARS-UCD1.2 chr7:90845941, BovineHD0700027239, BCHF5) and NFIA intron 4 (ARS-UCD1.2 chr3:84578325, BovineHD0300024307, BCHF31). However, genotyping tests involving tightly linked SNPs were also used on occasion. The assay steps consisted of sample processing (blood and ear notch), target amplification, and detection. Samples were processed and analyzed on the Solas 8 instrument per the manufacturer’s instructions. Briefly, samples were mixed with a proprietary lysis buffer and heated for three minutes at 95°C. Kit reagents were added to perform an initial PCR amplification with padlock probes. The machine assay conditions were SNP-specific and uploaded to the machine via the internet by the supplier: an initial heating phase at 93° for 3 min, followed by 30 cycles at 93° for 20 s and 57°C for 30 s. Detection of the ligated circular oligonucleotide products in the Solas 8 instrument was accomplished with fluorescently labeled probes and isothermal rolling circle amplification. The total reaction time was approximately two hours, and the genotype assignments were made in real time by the instrument, and the real-time graphical outputs were visually reviewed for quality assurance. Assay performance was validated with known genotypes derived from whole genome sequencing and the BovineHD BeadChip array from a beef cattle diversity panel (Heaton et al., 2016). Ambiguous genotype results from animals of unknown genotypes were resolved by repeat testing with a second assay designed for the opposite strand of ARRDC3 C182Y (ARS-UCD1.2 chr7:90845941, BovineHD0700027239, BCHF5r) and the linked NFIA intron 4 SNP (ARS-UCD1.2 chr3:84580655, BovineHD0300024308, BCHF32).

Genotype scoring, SNP filtering, and population substructure

Approximately 778,000 SNPs were scored for 204 animals in matched case-control pairs with a total genotyping rate of 0.99 (Files S1 and S2, Extended data). Subsequent filtering removed 216,000 SNPs with either low minor allele frequencies, missing genotypes, or deviation from Hardy-Weinberg expectations (173,000, 42,000, and 1,000, respectively). Approximately 560,000 remaining SNPs (Files S3 and S4, Extended data) were used for MDS analysis of pairwise IBS distances of the 204 animals from the 102 matched pairs, together with a 19-breed beef diversity panel. The feedlot animals from the 102 matched pairs were tightly clustered with purebred American Angus and Red Angus cattle (Figure 2A). A higher resolution plot with the 102 pairs of BCHF cases and matched controls showed there were seven outliers from one feedlot group, of which five were cases and two were controls (Figure 2B, lower right quadrant). However, many of the matched pairs had short IBS distances, indicating close genetic relationships between random penmates, as would be expected from animals originating from the same ranch. Based on the relatively low number of outliers that were BCHF cases, and the distances of their matched pairs, corrections for population stratification were not made in subsequent analyses.

Figure 2. MDS plot with pairwise IBS distances between breeds and paired individuals.

Panel A, genetic distances between BCHF case-control animals and diverse purebred bulls from 19 U.S. beef breeds. Panel B, genetic distances between matched pair from different feedlot groups and locations in Nebraska (NE) and Wyoming (WY). Lines connecting dots identify each matched pair with numbers corresponding to the pair ID (Table S1) and the clinical case.

Genome-wide McNemar association testing

The CMH test was first used to analyze the data set with PLINK software, since it was available and tests for SNP-disease association conditional on clustering of the individuals and is a generalization of McNemar’s test. The CMH test is expected to give the same results as McNemar’s test when each cluster contains a pair of individuals since their test statistics are identical when each stratum shows a pair. However, CMH output in PLINK did not provide easy access to the individual McNemar’s quadrants results for each SNP nor the ability to test 1-copy, 1 or 2-copy, and 2-copy risk allele models. Thus, we sought to develop a custom process on a programming and numeric computing platform as described in the Methods section. After producing the tables of McNemar test statistics (File S11, Extended data), 704, 668, 31, and 1 SNPs were further removed because they had no informative pairs, were unmapped, on the chromosome Y, or on the mitochondrial genome, respectively (File S12, Extended data). A Manhattan plot of the resulting SNPs revealed that the most significant genome-wide associations were clustered on chromosome 7 at the ARRDC3 gene (Figure 3). No other genomic region contained SNPs that met the Bonferroni threshold at the 0.05 genome-wide significance level (-log₁₀ = 7.05).

Figure 3. Manhattan and Q-Q plots of genome-wide association with 102 BCHF matched case-control pairs.

The mid-p-values from the McNemar’s association test were used with major (left) and minor alleles (right). In both panels, the inset shows the Q-Q plot of the distribution of the test statistics and the red horizontal line shows the Bonferroni correction at an α of 0.05.

While the Bonferroni correction protects against Type I errors (false positives), it is vulnerable to Type II errors (false negatives). To strike a balance between the number of true and false positives, q-values were calculated for the McNemar data set to identify regions surpassing a 5% FDR. The q-value threshold is the expected proportion of significant SNP associations that are false leads. There were 53 genome-wide autosomal SNPs passing the 5% FDR: 49 were associated with ARRDC3 on chromosome 7, five were associated with NFIA on chromosome 3, and one SNP each was associated with a region on chromosomes 4, 10, and 19 (Figure 4, File S13 Extended data). Since the three SNPs on chromosomes 4, 10, and 19 correspond to the expected proportion of false leads, these regions were dismissed for the purposes of the present study. The remaining genomic regions associated with ARRDC3 and NFIA genes were further characterized.

Figure 4. Genome-wide SNPs with a q-value threshold less than 0.05 FDR.

The q-values for SNPs were evaluated for the two-copy risk allele model and plotted by chromosome. Numbers in parentheses are linked SNPs with q-values less than 0.05. Red dots and parentheses represent multiple linked SNPs with equivalent q-values.

The best associated region spanned 52 kb and was roughly centered on the ARRDC3 gene. This region contained eight equivalent linked SNPs with -log₁₀ mid-p-values of 7.86 (Figure 5A) and contained a missense variant in ARRDC3 codon 182 in exon 4 encoding cysteine (C, tgt) or tyrosine (Y, tat). The homozygous major alleles for the eight best SNPs (and ARRDC3 Y182) were associated with increased BCHF risk (OR 8.4, CI₉₅, 3.3 to 21, Table 2). Conversely, one or two copies of the minor alleles for these eight SNPs (and ARRDC3 C182) were associated with reduced BCHF risk (OR = 0.12, File S12, Extended data). Heterozygosity (i.e., having one risk allele) was also associated with reduced BCHF risk (0.26). The second best associated region spanned a 2.3 kb region in NFIA intron 4, and there were no SNPs present on the BovineHD BeadChip array predicted to affect the NFIA coding sequence. The 2.3 kb region contained two equivalent linked SNPs with -log₁₀ mid-p-values of 6.60 (Figure 6). The homozygous major alleles for these two SNPs were associated with increased BCHF risk (OR 7.4, CI₉₅, 2.9 to 19, Table 2). Conversely, one or two copies of the minor alleles for these two SNPs were associated with reduced BCHF risk (OR = 0.14, File S12 Extended data). Heterozygosity was also associated with reduced BCHF risk (0.40). Thus, two copies of the major ARRDC3 and NFIA haplotype alleles were associated with BCHF risk.

Figure 5. Manhattan plot of ARRDC3 gene region.

The mid-p-values from the McNemar’s association test of the most frequent SNP alleles were used. Panel A, test requiring the presence of homozygous risk alleles. Panel B, test requiring the presence of either 1- or 2-copies of the risk alleles. Panel C, test requiring heterozygosity for the risk alleles. The red horizontal line shows the Bonferroni correction at the 0.05 significance level. The red numbers in parentheses are McNemar's informative pairs from quadrants b and c (i.e., “b/c“, Table 1). The ratio of informative pairs indicates the strength of the association, and whether the allele is positively or negatively associated with BCHF (i.e., b/c > 0 or < 0).

Figure 6. Manhattan plot of NFIA gene region.

The mid-p-values from the McNemar’s association test of the most frequent SNP alleles were used. Panel A, test requiring the presence of homozygous risk alleles. Panel B, test requiring the presence of either 1- or 2-copies of the risk alleles. Panel C, test requiring heterozygosity for the risk alleles. The red horizontal line shows the Bonferroni correction at the 0.05 significance level. The red numbers in parentheses are McNemar's informative pairs from quadrants b and c (i.e., “b/c“, Table 1). The ratio of informative pairs indicates the strength of the association, and whether the allele is positively or negatively associated with BCHF (i.e., b/c > 0 or < 0).

Conservation of ARRDC3 cysteine codon at position 182

The association of a missense variant in ARRDC3 with BCHF raised the possibility that the amino acid substitution could affect protein function and thereby increase the risk of disease in feedlot cattle. A comparison with ARRDC3 codons at position 182 in the jawed-vertebrate species (Gnathostomes) showed that cysteine residues were invariant throughout the 41 representative species examined (Figure 7). Examples of Gnathostome species with alternative residues at the equivalent position were not observed in any genus. The conservation of the C182 residue throughout the Gnathostomes is consistent with the hypothesis that the C182 residue is critical for normal ARRDC3 protein function and that homozygosity of the B. taurus-specific Y182 variant may have a deleterious effect.

Figure 7. Comparison of ARRDC3 residues near the C182Y position in the Gnathostomes.

TMRCA is the estimated time to the most recent common ancestor in millions of years (Hedges et al., 2015). The full length ARRDC3 protein is 414 residues in cattle and most of the Amniota species. The letters are IUPAC/IUBMB codes for amino acids. The dots are amino acid residues identical to those in cattle.

Predicting BCHF risk with ARRDC3 and NFIA SNP genotypes

A custom two-SNP test was designed on a commercially available platform to genotype cattle for the major ARRDC3 and NFIA risk alleles. The platform uses “padlock” oligonucleotide probes, combined with fluorescently labeled probes and isothermal rolling circle amplification technology, for the qualitative detection of alleles. Using reference DNA from animals with known genotypes, the call rate and accuracy typically exceeded 97% and 99%, respectively. This included repeated testing of approximately 10% of the samples due to ambiguous allele calls on the first pass. Using custom genotypes to resolve errors in the original BovineHD BeadChip array, the BCHF OR for calves with homozygous risk alleles at both genes increased to 28-fold more (b/c = 28/1, p-value = 1.1 × 10^-7, CI₉₅ = 4-206). Testing an independent set of 171 BCHF cases from the same feedlots showed approximately 29% were homozygous for both ARRDC3 and NFIA risk alleles, compared to 29% of the BCHF 102 cases used in the GWAS (Figure 8). Conversely, an independent set of feedlot cattle persistently infected with BVDV had about 2% homozygosity for both ARRDC3 and NFIA risk alleles, which was similar to the 3% observed with the 102 unaffected animals used in the GWAS. Assuming a 5% BCHF prevalence in the most severely affected pens of feedlot cattle, the adjusted negative and positive predictive values were 0.99 and 0.10, respectively, with a sensitivity and specificity of 0.91 and 0.58. Thus, genetic testing with two SNPs has significant power to identify which animals have the least genetic predisposition for developing BCHF in feedlots, while having less power to predict which animals will develop disease. ARRDC3 and NFIA genotype configurations and their corresponding expected disease and allele transmission risks, together with their allele frequency estimates by breed, are presented in Tables 3 and 4, respectively. Although imperfect, these genetic associations and DNA tests provide a starting point for reducing BCHF risk in the most affected herds.

Figure 8. Distribution of homozygous ARRDC3 and NFIA risk alleles in cattle cohorts.

Two-SNPs were used for identifying feedlot cattle with homozygous risk alleles (BovineHD0700027239 [ARRDC3] and BovineHD0300024308 [NFIA]). Bovine cohorts: USMARC beef cattle diversity panel v2.9 (19 breeds (Heaton et al., 2016), Kansas feedlot calves persistently infected with bovine viral diarrhea virus (BVDV PI (Workman et al., 2016), pen matched BCHF controls (this study), BCHF cases (this study), independent BCHF validation cases (this study). The BCHF validation cases were collected by pen riders experienced in BCHF diagnosis, however, they were not confirmed at necropsy by veterinarians and researchers. All 204 of the BCHF case-control cattle tested negative for BVDV, and thus were not persistently infected with BVDV.

Underlying and extended data

Figshare: Additional files for the article entitled, Association of ARRDC3 and NFIA variants with bovine congestive heart failure in feedlot cattle, https://doi.org/10.6084/m9.figshare.c.4331114.v2 (Heaton, 2022).

This project contains the following data files:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

Analysis code

MatLab analysis code available at GitHub: https://github.com/greg-harhay/McNemarsSNPAnalysis

Executable MatLab code available at Code Ocean: https://doi.org/10.24433/CO.0870729.v1

Data are available under the terms of the UnLicense.